Sofiane Bousabaa

M Arnaud

Derode Examinateur

M Philippe

Druault Examinateur

M Sandrine

Fauqueux Examinateur

M Emmanuel

Julliard Examinateur

M Sofiane

Khelladi Examinateur

M Régis

Marchiano Directeur De Thèse

M Daniel-Ciprian

Invité M Mincu

Ollivier François

M Invité

Radoi Emanuel

M Rapporteur

Valeau Vincent

Rapporteur

Sorbonne Universit

R Ésum

Arnaud Derode

Philippe Druault

Emmanuel Julliard

ACOUSTIC GREEN'S FUNCTION ESTIMATION USING NUMERICAL SIMULATIONS AND APPLICATION TO EXTERN AEROACOUSTIC BEAMFORMING

Acoustic imaging techniques aims at characterizing the different acoustic sources of noise on an aircraft (jet, fans, wing, landing gears, etc.) using microphone array measurements. Those source identification techniques, among whom the beamforming is the most popular, require the knowledge of the acoustic Green's function of the medium between the estimated sources and the microphones of the array. Unfortunately, those propagation functions are known only for cases of relatively simple complexity. In the presence of complex flows, cavities, reflecting surfaces like wings and tails, the Green's function is not known and the use of beamforming techniques based on approximate Green's function can lead to errors in the estimation of the location and amplitudes of the sources and can even lead to the apparition of spurious sources. The main aim of this thesis is to set up a numerical method for the estimation of the Green's function for aeroacoustic imaging applications. The method must have a minimal computational cost and provide a sufficiently accurate estimation to be used on realistic industrial configurations. The proposed methodology takes advantage of the sparsity of the Green's functions in the time-domain to minimize the required simulation time. The close relationship with the domain of system identification makes possible the use of a wide variety of sparsity-based regression algorithms like, among others, the stepwise regression, the lasso and the elastic net. First, the method is validated on complex 3D numerical test cases in the presence of flows and diffracting objects that are typical of those encountered in the industry. For the configurations involving a high number of focus points, the reverse-flow reciprocity decreases significantly the Green's function estimation problem. The methodology is finally applied on experimental data obtained on a high lift 2D wing placed in the ONERA CEPRA19 open section anechoic wind tunnel justifying the applicability of the method on realistic industrial configurations.

Introduction

The increase in the air transports during the last decades has made the aircraft noise reduction a necessity. During the phase of conception of aircrafts, pilots ensure that, during the take-off and landing phases, a minimum of noise is radiated in the direction of the ground. On the other hand, engineers and researchers work in creating quieter aircrafts and propellers [START_REF] Delrieux | From design to flight testing: overview of rotorcraft acoustic research at onera for industrial applications[END_REF]; [START_REF] Huber | Large-scale jet noise testing, reduction and methods validation "exejet": 1. project overview and focus on installation[END_REF]). It is thus important to give industrials a robust and precise methodology to enable the localization of acoustic sources on the wind tunnel model during the wind tunnel tests phase (cf. figure 1).

Acoustic source localization techniques aim to detect, localize and characterize acoustic sources on a given area of interest using an array of microphones that collects spatial and temporal samples of propagating wave fields. The set of acquired microphone signals are then processed in order to extract the relevant information. Several methods can then be used to recover sound sources: holography [START_REF] Ahmed | Holography and its application to acoustic imaging[END_REF]) beamforming [START_REF] Soderman | Directional microphone array for acoustic studies of wind tunnel models[END_REF]; [START_REF] Billingsley | The acoustic telescope[END_REF]; [START_REF] Elias | Experimental techniques for source location[END_REF]), and time reversal [START_REF] Tourin | Timedependent coherent backscattering of acoustic waves[END_REF]; [START_REF] Fink | Time-reversed acoustics[END_REF]; [START_REF] Deneuve | A coupled timereversal/complex differentiation method for aeroacoustic sensitivity analysis: towards a source detection procedure[END_REF]; [START_REF] Padois | Experimental localization of an acoustic sound source in a wind-tunnel flow by using a numerical time-reversal technique[END_REF]) among others. Beamforming methods require the knowledge of the acoustic impulse response of the medium or the Green's function (GF). However, those functions are known analytically on small number of configurations and under numerous approximations. In the presence of complex geometries and flows, the GF has no analytical expression. In practice, there are three possible ways for obtaining the GFs: analytical, experimental and using numerical simulation. Considering the analytical approach, [START_REF] Amiet | Refraction of sound by a shear layer[END_REF]; [START_REF] Allen | Image method for efficiently simulating smallroom acoustics[END_REF]; [START_REF] Bowman | Electromagnetic and acoustic scattering by simple shapes[END_REF]; [START_REF] Rienstra | An analytic green's function for a lined circular duct containing uniform mean flow[END_REF]; [START_REF] Sijtsma | Greens functions for in-duct beamforming applications[END_REF] have determined the GF for geometries up to annular ducts with a stepwise constant flow in the radial direction. However, in complex media with flows, scattering objects and reflecting surfaces, an analytic expression of the GF remains difficult to

obtain. An alternative is the experimental approach, [START_REF] Malbequi | Wind tunnel study of 3D sound propagation in presence of a hill and of a wind gradient[END_REF]; [START_REF] Koop | Investigation of the systematic phase mismatch in microphone-array analysis[END_REF] consists in placing a calibrated acoustic source for every focus point to determine the effects of the acoustic propagation. Two important disadvantages of this approach come from the difficulty to create an acoustic source that is isotropic and that has a low intrusivity. These two drawbacks can be addressed by using a laser-based sound source as envisaged recently by [START_REF] Rossignol | Validation of DLR's sound shielding prediction tool using a novel sound source[END_REF]; [START_REF] Rossignol | Analysis of the noise shielding characteristics of a NACA0012 2D wing[END_REF]. However, this solution is not an option in view of the wind tunnel operating costs and filled order books. The third approach is to determine the GF by using numerical simulations and is the one we are interested in in this study. It has the advantage to be less limiting in terms of complexity of both the geometry and the flow. Moreover, the progress in Computational AeroAcoustics (CAA) enables to consider realistic configurations with low numerical dissipation [START_REF] Tam | Computational aeroacoustics -issues and methods[END_REF][START_REF] Lele | Computational aeroacoustics -a review[END_REF]. However, the extra precision go hand in hand with an increase of the computational cost. [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationary computational fluid dynamics[END_REF] proposed an approach using excited Computational Fluid Dynamics simulations with system identification (CFD/SI) to extract acoustic transfer functions of systems containing two or three inputs [START_REF] Föller | Aeroacoustic characterization of T-junctions based on large eddy simulation and system identification[END_REF]. The system is modeled as a Linear Time Invariant (LTI) system and system identification [START_REF] Ljung | System identification[END_REF] techniques are applied to time series data to extract the unit impulse responses of the multi-port.

To do so, they compute the auto-and cross-correlations of the input signals and the impulse response of the system is obtained by the inversion of a Wiener filter equation.

This approach enables an estimation of the impulse response for all the input/output pairs from one single simulation. However, this approach requires the knowledge of the reverberation time and the computational cost is increased when this time is big or overestimated.

The aim of this thesis is to develop a methodology for the estimation of the GFs based on numerical simulations for extern aeroacoustic imaging applications. The method must have a minimal computational cost and provide Green's functions that are sufficiently accurate to be used on realistic industrial configurations.

It represents a continuation of the one of [START_REF] Pene | Antennerie numérique pour la caractérisation de sources aéroacoustiques en milieu complexe[END_REF]. He arrived to an ill-posed problem that he decided to regularize by using a constraint on both the l 1 and l 2 norm of the solution. The resulting methodology allowed the estimation of the GF on a number of 2D test cases like the diffraction by a disc or by a half-plane. However, it presented some handicaps, in particular a high post-processing cost and the a priori knowledge of some informations on the GFs for the calibration of the stopping criteria.

The manuscript is organized as follows. The chapter 1 is dedicated to the definition of the acoustic Green's functions and their role in acoustic imaging. First, the mathematical background of the acoustic GF is exposed. We show that it is the impulse response of a linear system representing the propagation through the medium. Acoustic imaging techniques are then presented and the role played by the GF and its influence on the quality of acoustic maps will be highlighted. The various approaches for the obtention of the GFs are then described. The advantages and drawbacks of analytical, experimental, and numerical approaches are discussed and we show that, to reduce the computational cost, one can be tempted to estimate all the required Green's functions from a single CAA simulation.

The proof of the ill-posed nature of the GF estimation problem is given in the chapter 2. In order to regularize it, it is necessary to add a priori information on the solution to the model. In the case of external aeroacoustic imaging, in the absence of resonance or reverberation, the acoustic GFs can often be considered sparse. We show that the adding of a constraint on the l 0 , l 1 , or l 2 -norm of the solution enables to regularize the GFs estimation problem. Several algorithms for taking into account the sparsity are adapted and compared. Those iterative algorithms require a stopping criterion and the cross-validation is presented for this purpose.

In the chapter 3, a reflexion on the choice of the signals is led. An iterative algorithm is proposed for the generation of multisine signals [START_REF] Schoukens | Survey of excitation signals for FFT based signal analyzers[END_REF] with low crosscorrelation. It is based on the crest-factor minimization algorithm found in [START_REF] Guillaume | Crest-factor minimization using nonlinear Chebyshev approximation methods[END_REF], approximating the Chebyshev norm with l p -norms with increasing values of p. We also propose an adaptation of the algorithm to the case of finite impulse response (FIR) systems whose estimation requires the cross-correlation to be low on a limited time interval. The contributions of the cross-correlation minimization to the GF estimation are finally discussed.

In the chapter 4, two test cases representative for extern shapes for fuselage and wings are considered: the diffraction by a rigid sphere and the diffraction by a NACA0012

wing profile. For the first test case, the pressure field resulting from the diffraction of a monopole source by a rigid sphere can be expressed analytically. In order to prove the advantages of CAA in terms of GF estimation, we will compare CAA estimated GFs with both experimental and analytical ones. For the second test case, the solution is not known analytically and the GFs obtained using the proposed method will be compared only with experimental ones. Beamforming results will also be presented proving that the method can be applied to real industrial configurations with several thousands of focus points.

Finally, the method is applied in experimental microphone array data measurements in the chapter 5. Data were collected in the ONERA CEPRA19 anechoic wind tunnel [START_REF] Piccin | Cepra19: The ONERA large anechoic facility: A major tool for aeroacousitic measurements[END_REF] and were part of the European Commission (EC) Technologies to IMProve Airframe Noise (TIMPAN) project [START_REF] Manoha | Slat noise reduction using meshes: acoustic tests in CEPRA19 and numerical simulation of the 3D flow (EC project TIMPAN)[END_REF][START_REF] Perrin Decroux | TIMPAN -Technologies to IMProve Aiframe Noise[END_REF]. The mean flow is obtained using CFD simulation solving the RANS equations and is then used during the CAA simulation in subsection 5.2.2. Beamforming maps using CAA estimated Green's functions are compared with ones obtained using classical shear layer correction method.

Chapter 1

Methods for the obtention of the acoustic Green's function

This chapter is dedicated to the definition of the acoustic Green's functions and their role in acoustic imaging. First, in section 1.1, the mathematical background of the acoustic GF is exposed. We show that it is the impulse response of a linear system representing the propagation through the medium. Acoustic imaging techniques are then presented in the section 1.2. The role played by the GF is exposed and its influence on the quality of acoustic maps. Finally, the various approaches to obtain the GFs are considered in section 1.3. Advantages and drawbacks of analytical, experimental, and numerical methods are discussed.

1.1 The acoustic Green's function

Acoustic model

Euler equations are a set of nonlinear partial derivative equations (PDE) that describe the flow of a fluid or gas. They correspond to a particular case of the Navier-Stokes equations in the absence of viscosity and thermal conduction. These equations are a good model for the propagation of the sound and can be solved numerically by the use of a computer. Details on Euler equations and their numerical implementation are given in [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]. However, aeroacoustics simulations using Euler equation is still a complex problem because it involves different scales. The propagation of acoustic waves through a moving medium can be modeled by the linearized isentropic Euler equations [START_REF] Bailly | Numerical solution of acoustic propagation problems using linearized euler's equations[END_REF]. Euler equations are linearized around a given stationary mean flow. This mean flow can be obtained analytically or by solving a simplified problem not involving acoustics (Reynolds-averaged Navier-Stokes equations, potential

Chapter 1 Methods for the obtention of the acoustic Green's function flow equations, . . . ). Linearized Euler system of PDE reads as:

       ∂p ∂t + u 0 • ∇p + u • ∇p 0 + γp 0 ∇ • u + γp∇ • u 0 = -c 2 0 Θ ∂u ∂t + u 0 • ∇u + u • ∇u 0 + 1 ρ 0 ∇p -p ρ 2 0 c 2 0 ∇p 0 = 0 p = c 2 0 ρ and c 2 0 = γp 0 ρ 0 (1.1)
The density ρ, the velocity u and the pressure p designate small perturbations superimposed on a mean steady flow of density ρ 0 , velocity u 0 and pressure p 0 . The coefficient γ designates the ratio of specific heats, and is taken as γ = 1.4 for air. The term Θ at the right hand side

-c 2 0 Θ(x, t) = -c 2 0 ∂Θ ∂t
is used to excite the medium and represents an injection of mass.

The GF g (y,τ ) (x, t) is the pressure field resulting from a pulse both in spatial (at position y) and temporal domain (at time τ ). In other words, it is the solution of (1.1) for pressure when the forcing term is Θ(x, t) = δ(xy)δ(tτ ). The Dirac delta function is the identity element for the convolution and we have the property

Θ(y, τ ) = τ -∞ x∈V Θ(x, t)δ(x -y)δ(t -τ ) d 3 x dt,
and the superposition principle enables to write p(x, t) = t -∞ y∈V Θ(y, τ )g (y,τ ) (x, t) d 3 y dτ.

(1.2)

If the mean flow is time-independent, the system is said to be linear time invariant (LTI): if the causes are delayed by an interval of time τ , the consequences (pressure and flow perturbations) will also suffer from the same delay. Mathematically, it means that g (y,τ ) (x, t) = g (y,0) (x, tτ ) and the integral in time in the equation (1.2) can be seen as a convolution product p(x, t) = s i (•) * g (y i ,0) (x m , •) (t).

(1.5)

This enables to separate the sound generation of the source from the propagation in the medium.

Signal p m (t) is the measure of the pressure fluctuations at the position x m of microphone m and the time function g i,m (t) = g (y i ,0) (x m , t) will be referred to as the GF between source i and microphone m. It represents the acoustic impulse response of the medium between those two points

p m (t) = N S i=1
(g i,m * s i ) (t) ∀m, t.

(1.6)

For each microphone m, this is the equation of an LTI Multiple Input Single Output (MISO) system as shown in the figure 1.1. This new formulation exhibit the time GFs as the impulse responses of the system and system identification (SI) [START_REF] Ljung | System identification[END_REF] techniques can be used for their determination.

Free-field solution

In a medium at rest (u 0 = 0, ∇p 0 = 0), the linearized Euler system of equations can be decoupled. After performing the derivation of the first equation in time, it is now possible to obtain a single PDE for the pressure fluctuations. To determine the freefield GF, a dirac delta function in the 3D space is placed in the right hand side. The invariance with respect to translations enable to place the source at the origin of the space:

       ∆g (0,0) - 1 c 2 0 ∂ 2 g (0,0) ∂t 2 = δ(x)δ(t) c 2 0 = γp 0 ρ 0 .
(1.7)

The Green's function G (0,0) is solution of:

       1 r 2 ∂ ∂r r 2 ∂g (0,0) ∂r - 1 c 2 0 ∂ 2 g (0,0) ∂t 2 = δ(x)δ(t) ∀(x, t) ∈ Ω × R + boundary conditions (1.8)
By performing a Fourier transform:

F (x, ω) = t∈R f (x, t)e -iωt dt (1.9)
Problem (1.8) becomes:

       1 r 2 ∂ ∂r r 2 ∂G (0,0) ∂r + k 2 G (0,0) (r, ω) = δ(x) ∀(x, t) ∈ Ω × R + boundary conditions (1.10)
For r > 0, the term in the right hand side of the first equation in (1.10) is zero. The solution space is the vector space of dimension 2 spanned by the functions e ikr /r and e -ikr /r. From those two functions, only the last one satisfies the Sommerfeld [START_REF] Orlanski | A simple boundary condition for unbounded hyperbolic flows[END_REF] radiation condition and thus G (0,0) (r, ω) = A e -ikr r .

(1.11)

To determine the constant A, we replace the expression in the equation (1.10), we multiply both sides by a test function φ and integrate over space (1.12)

It leads to -4πAφ(0) = φ(0) =⇒ A = - 1 4π . (1.13)
With the chosen convention, the frequency and time domain Green's functions read as G (y,0) (x, ω) = -e -i2πf x-y /c 0 4π xy and g (y,0) (x, t) = -δ (txy /c 0 ) 4π xy . (1.14)

It is expressed in kg/s 2 = (kg * m/s 2 )/m = N/m = P a * m. Since the 3D free-field GF decreases as 1/r, this results in the intensity of the acoustic wave decreasing as 1/r 2 and reflects the conservation of energy for spherical waves.

Role in acoustic imaging

The increased restrictions on aircraft noise during takeoff and approach have led to an increased interest in the localization of noise sources and quantification of their levels.

Hydrophone arrays [START_REF] Akyildiz | Underwater acoustic sensor networks: research challenges[END_REF] have been used for many years in sonar applications to detect ships at large distances. Because of the distance, acoustic waves arriving at the microphones may be considered as plane waves. In the following, L is the array overall typical size and λ the acoustic wavelength. Two cases can occur depending on the distance between the sources and the microphone array D. On one hand, we say that the sources are in the Fraunhofer region of the array if D L 2 /λ. In that case, only the directions of sources are looked for instead of their spatial location (See figure 1.2-(a)). One the other hand, the sources are said to be in the Fresnel region of the array if D ≤ L 2 /λ. The spherical characteristics of the wavefronts are not lost and the location of the acoustic source can be retrieved at the center of the spherical waves (Figure 1.2-(b)). Methods can be used in wind tunnel measurements or in real aircraft noise identification using large aperture phased arrays [START_REF] Brusniak | Acoustic imaging of aircraft noise sources using large aperture phased arrays[END_REF]. It is then possible to test different strategies for aircraft noise reduction [START_REF] Piet | Flight test investigation of add-on treatments to reduce aircraft airframe noise[END_REF].

Beamforming in the Fraunhofer region

The figure 1.2-(a) represents the case of acoustic imaging in the Fraunhofer region. The acoustic source is located far from the array in the direction u s . Resulting wavefronts are orthogonal to this vector. The signal registered at microphone m is expressed as

p m (t) = p 0 t + u s • x m c 0 ,
with p 0 (t) the pressure signal at the reference microphone located at the center of the array and at the origin of the spatial system and c 0 is the speed of sound. Thus, having the pressure fluctuations at microphones, the direction of the source can be retrieved.

To do so, the sum of these signals with compensation delays is computed

P ABF (u, t) = N M m=1 w m p m t - u • x m c 0 .
The angular beamforming operator P ABF (u, t) is used to detect if wavefronts are coming from the direction u. The previous formulation is often dropped in favor of the frequency domain formulation .15) with p m (f ) the complex amplitude of pressure fluctuations at microphone m. Even if this expression was derived for a unique source, there may be several sources in presence. The quantity P ABF (u, f ) is computed for all the directions u. High values of this function are obtained in directions were sources are more likely to be.

P ABF (u, f ) = N M m=1 w m p m (f ) exp -2iπf u • x m c 0 , ( 1
One way to study the behavior of an acoustic array is to consider the case of a line array and a point monopole source. The microphones are taken equally spaced by a distance δ in the X-axis, x m = mδu x . If one source is emitting signal from the direction u s , the pressure at the microphones is

p m (f ) ∝ exp 2iπf mδ u s • u x c 0 .
If all the receiver weights are chosen unitary, w m = 1 for all m, the beamformer output 1.15 can be computed analytically and we obtain

|P ABF (u, f )| ∝ sin πf N M δ(u-us)•ux c 0 N M sin πf δ(u-us)•ux c 0 .
This function is plotted in the figure 1.3-(a) at the frequency 1 kHz, 2 kHz and 4 kHz.

The line array contains N M = 25 microphones and the distance between two neighbor microphones is δ = 1 cm. The medium is air at rest and at ambient temperature and the speed of sound is c 0 = 340 m/s. The spatial resolution, or selectivity, of the beamformer is the capacity to identify two sources close to one another. It is related to the width of the main lobe and the amplitude of the side lobes that must be as small as possible.

This width is inverse proportional to the considered frequency it can be seen in the figure 1.3-(a). The beamforming output at 1 kHz is shown in black color, the width of the main lobe is divided by 2 at 2 kHz and by a factor 4 at 4 kHz. Spatial aliasing can occur if the distance between two neighbor microphones is too big.

Very high side lobes appear and can be an obstacle to the determination of the true direction of the source. This effect is visible in the figure 1.3-(a) at 4 kHz. The Shannon criterion provides conditions to avoid spatial aliasing. This criterion is respected if the distance between two neighbor microphones is smaller than half the wavelength.

The choice of the receiver weights w m has a significant effect on the shape of the ar-

Beamforming in the Fresnel region

In aeroacoustic measurements [START_REF] Mueller | Aeroacoustic measurements[END_REF], the source is often located at close distance to the array in the Fresnel region of the array D ≤ L 2 /λ. This case is also referred to as focused beamforming and is represented in the figure 1. 2-(b) . Waves arriving at microphones cannot be considered as plane anymore. Strictly speaking, the previous results on lobes characteristics and spatial aliasing are not exact. However, it appears in practice that the same behavior is obtained with focused processing.

In a medium at rest and without diffracting objects, the signal registered at microphone m is expressed as

p m (t) = 1 4πr ms s t - r ms c 0 ,
with s(t) is the source signal and r ms the distance between the source and the microphone. This expression was derived from the case of one monopole source, however the position of the source is not known a priori. The application of the beamforming method consists in several steps. First, the zone containing the sources is guessed and discretized into N S focus points. Then, for each focus point i, a beamformer output P DAS (x s , t) can be formed by averaging weighted and delayed versions of the receiver signals

P DAS (x i , t) = N M m=1 w m p m t + x m -x i c 0 . (1.16)
This procedure is called delay-and-sum (DAS) and it is a time-domain beamforming operator. The delays compensate the propagation times from sources to microphones and high output values are obtained at the positions of the sources. It was shown by [START_REF] Fischer | Beamforming of aeroacoustic sources in the time domain[END_REF] that this technique allows assessing the intermittent nature of aeroacoustic sources. Details on this operator and its application for discrete time beamforming can be found in the book of [START_REF] Dudgeon | Multidimensional Digital Signal Processing Prentice-Hall Signal Processing Series[END_REF].

The DAS beamforming is limited when the domain contains reflecting surfaces and complex flows. Alternative time-domain imaging methods were proposed for studying aeroacoustic sources, relying on the time-reversal principle [START_REF] Fink | Time-reversed acoustics[END_REF][START_REF] Deneuve | A coupled timereversal/complex differentiation method for aeroacoustic sensitivity analysis: towards a source detection procedure[END_REF][START_REF] Padois | Experimental localization of an acoustic sound source in a wind-tunnel flow by using a numerical time-reversal technique[END_REF][START_REF] Mimani | Multiple line arrays for the characterization of aeroacoustic sources using a time-reversal method[END_REF][START_REF] Rakotoarisoa | Time-domain delay-and-sum beamforming for time-reversal detection of intermittent acoustic sources in flows[END_REF]. It is a promising method operating in the time-domain used for solving inverse problems of sound source localization. To do so, the recorded acoustic pressure timehistory is reversed in time and is enforced at the microphone location. Acoustic waves are then back-propagated after the mean flow direction was numerically reversed. Authors have proposed alternatives for more complex configurations for instance with the rotating sources. To image rotating fan blade noise, [START_REF] Dougherty | Virtual rotating microphone imaging of broadband fan noise[END_REF] applied timedomain beamforming using a beamforming grid that rotates with the fan.

While time-domain is a robust method operating in the time domain for solving inverse problems of sound source localization, most of the experimental studies in wind tunnel tests are still done in the frequency domain. The various sources often exhibit a different spectral behavior and their location depend on the frequency of interest. In the Fourier domain, relation (1.6) reads as:

P m (f ) = N S i=1 G i,m (f )S i (f ) ∀m, f.
(1.17)

From this relation, several estimators can be derived for the estimation of the power of the sources. Most of them are based on the microphone cross-correlation matrix

Γ(f ) = Γ m,m (f ) 1≤m,m ≤N M with Γ m,m (f ) = P (k) m (f )P (k) m (f ) * .
The cross-spectra Γ m,m (f ) of the microphone signals are computed with the classical periodogram method [START_REF] Welch | The use of fast fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms[END_REF]. The recorded time signals are divided into N K blocks of data. The cross-spectrum is computed from the Fourier transforms of the microphone signals over each block k. The ensemble of cross-spectra are then averaged operation over the data blocks and overline denotes the statistical mean. In the periodogram method, there is a trade-off between the number of blocks and the frequency resolution. Cutting time signals into an high number of blocks will ensure a better statistical convergence but will result in larger frequency bins.

The conventional beamforming (CBF) operator is

P CBF (x i , f ) = m,m G * i,m (f )Γ m,m (f )G i,m (f ) m,m G i,m (f ) 2 G i,m (f ) 2 . (1.18)
This formula is widely used and enables to obtain acoustic maps that bring both source location and power. The steering vector

G i = [G i,1 (f ), • • • , G i,N M (f )] is the concatena-
tion of all the GFs for a focus point i. With this, the CBF operator takes the form

P CBF (x i , f ) = G + i (f )Γ(f )G i (f ) G i (f ) 4 , (1.19)
where superscript + is the hermitian operator. In his book, [START_REF] Mueller | Aeroacoustic measurements[END_REF] gives more details on the origins of this expression. If a unique source is present at focus point j and has an amplitude |S j (f )| at frequency f , the terms of the cross-correlation matrix will be Γ(f ) = G + j (f )G j (f ) and the CBF operator will give exactly the source power α j (f ) = |S j (f )| 2 at the source location x j . As it will be seen in section 1.2.3 and later in 4.2, this operator has a tendency to prioritize hidden focus points for which the GF amplitude |G i,m (f )| takes small values. The results are often plotted in dB with the formula

P dB CBF (x i , f ) = 20 log 10 P CBF (x i , f ) p 0 ,
with p 0 is a reference pressure solution that can be chosen as p 0 = 2 × 10 -5 Pa. It is near the absolute threshold for a normal human listener for a sound frequency of 1000 Hz.

In the following, CBF results will be plotted in dB. Sometimes, acoustic noise sources can exhibit broadband spectrum characteristics. For this type of sources, the considered frequency range is often discretized and acoustic maps are computed on all the resulting frequencies. The contributions of the various frequencies are then taken into account by summing the resulting acoustic maps.

An alternative method for source location is the coherence with the monopole model (CMM) [START_REF] Elias | Experimental techniques for source location[END_REF] and reads as

P CM M (x i , f ) = G + i (f )Γ(f )G i (f ) G i (f ) 2 Γ(f ) , (1.20)
where the euclidian norm on square matrices is given by M 2 = trace M + M . This operator gives values between 0 and 1. It is very precise in the localization of acoustic sources but doesn't give access to the amplitudes of the sources. The sources are more likely to be where this criteria takes high values and, if only one source is emitting in the absence of noise, the amplitude at the position of the source is equal to one.

! "#$ %&'#(% !"#$ !%#$ !&#$ !'#$ !(#$ #$ !($ !"($ !%($ !&($ !'($ (
a) CBF using free-field GFs.

! "## $ !"#$ !"%$ !"&$ !"'$ !$ ($ !")$ !"*$ !"+$ !",$ !"($ (

Influence of the accuracy of the Green's functions on acoustic maps

Beamforming maps in a medium at rest and using the free-field GF is presented in the figure 1.4. The 51 microphones are placed in the segment -0.5 m ≤ x m ≤ 0.5 m and y m = -0.5 m every 0.02 m. The source is placed at x s = (0, 0.5 m, 0) and its amplitude is 0 dB at the considered frequency f = 5 kHz. The CBF map is shown in the figure 1.4-(a). One main lobe is obtained at the position of the source surrounded by numerous secondary lobes. The spatial expanse of the main lobe is inverse proportional to the frequency and directly related to the resolution. The microphone array thus has a higher sensitivity in the angular than in the depth direction.

Acoustic waves are affected by the presence of a mean flow as it is shown in the figure 1.5.

The same geometry is used however, this time, a mean flow v = v u x with v = 100 m/s is present in the region y ≥ 0. The use of the free-field GF is not adapted as it can be seen in the figure 1.5-(a) and (b). It leads to an error in the localization of the source that is detected downstream its true position. [START_REF] Amiet | Refraction of sound by a shear layer[END_REF] proposed an analytical correction to the free field GF to take into account the refraction by a thin shear-layer.

The beamforming maps based on the resulting GF is shown in the figures 1.5-(c) and (d).

The use of the correct GF enables to retrieve the position and amplitude of the source.

The presence of reflecting/diffracting surfaces also have an influence on the propagation of sound. In the figure 1.6, beamforming is performed in the presence of a rigid sphere. The sensor array is composed of 41 microphones placed in the segment -0.4 m ≤ x m ≤ 0.4 m and y m = -0.4 m every 0.02 m. The 0.2 m-radius sphere is placed at the origin and the source position is now x s = (0.2 m, 0.3 m, 0) and the medium is at rest. Again, the free-field GF appears to be unadapted and acoustic diffracted by the sphere gives to the array the impression that the source is located at the right edge

! "#$ %&'#(% !"#$ !%#$ !&#$ !'#$ !(#$ #$ !($ !"($ !%($ !&($ !'($ (
a) CBF using free-field GFs.

! "## $ !"#$ !"%$ !"&$ !"'$ !$ ($ !")$ !"*$ !"+$ !",$
!"($ (b) CMM using free-field GFs.

! "#$ %&'#(% The source is placed at x s = (0, 0.5 m, 0) and its amplitude is 0 dB at the considered frequency f = 5 kHz. The medium is at rest in the bottom region (y<0) and the flow velocity is 100 m/s in the top region (y≥0). Acoustic maps obtained using CBF and CMM are shown using the free-field GF and a more precise GF that take into account the presence of the shear layer.

!"#$ !%#$ !&#$ !'#$ !(#$ #$ !($ !"($ !%($ !&($ !'($ (c) CBF using analytical GFs. ! "## $ !"#$ !"%$ !"&$ !"'$ !$ ($ !")$ !"*$ !"+$ !",$ !"($ (
Chapter 1 Methods for the obtention of the acoustic Green's function

! "#$ %&'#(% !" #$!" #%!" #&!" #'!" #(!"
(a) CBF using free-field GFs.

! "## $ !"#$ !"%$ !"&$ !"'$ !"($ !")$ !"*$ !$ !"+$ !",$ (
b) CMM using free-field GFs.

! "#$ %&'#(%

!" #$!" #%!" #&!" #'!" #(!"
(c) CBF using analytical GFs. of the sphere 1.6-(b). For this geometry, an analytic expression of the GF is available and is given in the book of [START_REF] Bowman | Electromagnetic and acoustic scattering by simple shapes[END_REF]. Acoustic CBF maps based on the analytic GF is shown in the figure 1.6-(c). The position of the source is still difficult to extract because the CBF operator prioritize focus points located in hidden region, that is to say, those located behind the sphere. The CMM operator gives a better estimation of the position of the source (see figure 1.6-(d)).

! "## $ !"#$ !"%$ !"&$ !"'$ !"($ !")$ !"*$ !$ !"+$ !",$ (
Approaches have been proposed to improve the accuracy of the beamforming results.

Conventional beamforming methods have proven their efficiency in source localization, but they are limited in the estimation of the power levels of extended sources. [START_REF] Blacodon | Level estimation of extended acoustic sources using a parametric method[END_REF] proposed an approach based on a spectral estimation method (SEM) in order to determine the actual sound-pressure levels of the acoustic sources found by the localization method. The method is successfully applied to compute the power levels of the main sources on a 1/11th-scaled model of the Airbus A320/A321 tested in CEPRA 19 wind tunnel. To improve medical ultrasound imaging results, [START_REF] Synnevag | Benefits of minimum-variance beamforming in medical ultrasound imaging[END_REF] proposed to find the weights in (1.16) by minimizing the variance of the beamformer output. However, all of those methods will benefit from a more accurate model of propagation of the sound for configurations involving complex flow and geometry, e.g., [START_REF] Fleury | Large-scale jet noise testing, reduction and methods validation "exejet": 5. analysis of jet-airfoil interaction noise by microphone array techniques[END_REF] analyzed the jet-airfoil interaction noise by microphone array techniques.

Deconvolution methods

We have seen in the previous section that the acoustic maps obtained by the previously described methods are blurred and present numerous sidelobes. In deconvolution problems, those maps are often referred to as "dirty" maps. They were first introduced in astrophysics where images obtained by telescopes are blurred due to diffraction by the diaphragm. Deconvolution algorithms use the known beam pattern, or point spread function (PSF), to distinguish between real structure and sidelobe disturbances. The PSF represents the beamforming map obtained in the presence of one monopole at a focus point j. This problem can be stated as a linear system Hα = b.

(1.21)

The vector of unknowns α, of size N S , contains the power of the sources. The N S × N S matrix H is the concatenation of all the PSF, the vector b contains the CBF values

(1.18) for each of the N S focus points:

H i,j = G + i (f )G j (f ) 2 G i (f ) 2 G j (f ) 2 and b i = G + i (f )Γ(f )G i (f ) G i (f ) 4 . (1.22)
It can be shown that, behind (1.21), lies the hypothesis of the sources being distributions of statistically independent noise radiators (incoherent sources).

Several algorithm have been proposed to estimate the power level α radiated by the sources. The CLEAN algorithm starts from the dirty map b and iteratively substract PSF corresponding to the focus point i 0 where b i 0 is maximum [START_REF] Högbom | Aperture synthesis with a non-regular distribution of interferometer baselines[END_REF]. This algorithm suppose some regularity in the PSF in particular that the maximum amplitude corresponds to the position of the source. This is not the case for the figure 1.6

where we saw that the maximum CBF amplitude is not at the position of the source. [START_REF] Sijtsma | Clean based on spatial source coherence[END_REF] proposed CLEAN based on spatial source coherence (CLEAN-SC). It performs a decomposition of the CSM into coherent components and thus, unlike other deconvolution methods, it does not use the PSF.

Chapter 1 Methods for the obtention of the acoustic Green's function [START_REF] Brooks | A deconvolution approach for the mapping of acoustic sources (DAMAS) determined from phased microphone arrays[END_REF] proposed the deconvolution approach for the mapping of acoustic sources (DAMAS). The DAMAS method naturally apply physically-necessary positivity constraint on the vector of solutions α that also make the problem more deterministic. The original DAMAS algorithm suffers from slow convergence but it is still widely used in the aeroacoustic community and in particular in wind tunnel tests. Several versions were then proposed among whom DAMAS2 and DAMAS3 [START_REF] Dougherty | Extensions of damas and benefits and limitations of deconvolution in beamforming[END_REF] that restrict the point spread function to a translationally-invariant, convolutional, form to improve the calculation time. DAMAS-C was proposed by [START_REF] Brooks | Extension of damas phased array processing for spatial coherence determination (DAMAS-C)[END_REF] for sources with spatial coherence. [START_REF] Fleury | Determination of acoustic directivity from microphone array measurements using correlated monopoles[END_REF] proposed two extensions to reduce the heavy computational cost of the DAMAS-C and DAMAS2-3 methods. In SC-DAMAS, [START_REF] Yardibi | Sparsity constrained deconvolution approaches for acoustic source mapping[END_REF] added a sparsity constraint to solve problem (1.21)

for a sparse distribution of the sources (when vector α has a small number of non-zero components).

Finally, the deconvolution problem (1.21) can also be seen as a minimization problem.

If the spatial distribution of acoustic sources is sparse, the problem can be formulated

as minimize { α 1 } subject to Hα = b. (1.23)
Several linear regression algorithms like lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], Least Angle Regression (LARS) [START_REF] Efron | Least angle regression[END_REF] or the elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] can then be used to solve this problem. The l 0 -norm can also be preferred to the l 1 -norm in the minimization but it will result in an acoustic map that is more sparse. [START_REF] Padois | Orthogonal matching pursuit applied to the deconvolution approach for the mapping of acoustic sources inverse problem[END_REF] proposed recently to apply the Orthogonal Matching Pursuit algorithm, that is based on l 0 -regularization, to the deconvolution approach for the mapping of acoustic sources.

Other regularization can also be envisaged, for instance, [START_REF] Dougherty | Extensions of damas and benefits and limitations of deconvolution in beamforming[END_REF] proposed a regularization based on Wiener filter.

1.3 The different approaches to obtain the Green's functions

In the previous section 1.2, we saw that the accuracy of the GF plays an important role in acoustic imaging. Beamforming based on approximative GF leads to errors in the determination of the position and amplitude of acoustic sources. In this section, we discuss the different analytical, experimental and numerical approaches to obtain the GFs.

Analytical methods

In most of the cases, conventional beamforming methods use the free-field analytical GFs to describe the acoustic propagation between noise sources and microphones. To do so, microphones must be placed in a way that the propagation of waves is not affected too much by the presence of the flow and obstacles. In the study of the noise inside a turbofan engine, [START_REF] Sijtsma | Using phased array beamforming to locate broadband noise sources inside a turbofan engine[END_REF] used a circular microphone array in the engine intake duct wall and the free-field GF is used for the definition of the steering vectors. Several methods have been proposed to deal with more complex situations. [START_REF] Amiet | Refraction of sound by a shear layer[END_REF] proposed an analytical correction to the free field GF to take into account the refraction by a thin shear-layer. This correction is widely used in the adjustments of the refraction of a cylindrical shear layer in open-section wind tunnel tests [START_REF] Elias | Correction des effets de réfraction dans la soufflerie sourde CEPRA19[END_REF][START_REF] Fleury | Analysis of jet-airfoil interaction noise sources by using a microphone array technique[END_REF]. For small-room acoustics, [START_REF] Allen | Image method for efficiently simulating smallroom acoustics[END_REF] investigated the use of image source techniques for simulating the impulse response between two points in a small rectangular room. This correction can also be applied to rectangular closed-section wind tunnel experiments.

In the 20th century, the electromagnetic waves by simple shapes interested researchers for the purpose of radar applications. In his book, [START_REF] Bowman | Electromagnetic and acoustic scattering by simple shapes[END_REF] gathers the analytical expression of the diffraction by several geometric shapes. Green's functions in the presence of finite or infinite body in two or three dimensions (disc, sphere, cylinder, cone, etc.). Unfortunately, those shapes appear inadequate to describe the complex shapes encountered in the industry. Another limitation of the analytical methods is the consideration of the influence of the mean flow on acoustics. [START_REF] Taylor | A transformation of the acoustic equation with implications for windtunnel and low-speed flight tests[END_REF] Recently, [START_REF] Rienstra | An analytic green's function for a lined circular duct containing uniform mean flow[END_REF] proposed an analytical GF in annular ducts has been calculated for a uniform flow. [START_REF] Sijtsma | Greens functions for in-duct beamforming applications[END_REF] extended it to the case of a stepwise constant flow in the radial direction. However, in the presence of scattering objects and reflecting surfaces, the GF is difficult to obtain especially in the presence of complex flows. Moreover, the analytical expression of the GF is also of increasing complexity.

The analytical expression of the GFs have become more and more complex to compute reducing the initial interest of this approach. Moreover, a lot of approximations are still to be made because of the complexity of cases encountered in the industry. The use of beamforming methods based on approximate GFs may result in biased sources positions and power estimates or even lead to spurious sources.

!"#$%&'%()*+$$+,*

p 1 p N M • • • • • • p m !)

Experimental methods

The analytical methods allow to determine the Green's function in relatively simple situations. However, for a number of cases of application encountered in aeroacoustics, there is no analytical solution. In the case of fan noise imaging on an installed turbofan, shown in the figure 1.7, one may be tempted to install microphones outside the fan cowl to reduce intrusivity. The propagation of acoustic waves is affected by the fan cowl and the flow making it difficult to localize acoustic sources. for an open wind tunnel using an in-flow calibration source. This source employs a ribbon loudspeaker which can provide sufficiently high sound pressure levels in a broad frequency range and its properties are known. Effects of the flow can be estimated and a correction is proposed to improve the quality of source maps and the estimation of source levels. While this study enables to test correction models in open wind tunnel, it is not adapted to GFs estimation because of the high intrusivity of the acoustic source and the high financial cost.

An alternative to reduce source intrusivity is to use a laser-based sound source. By focusing a high energy laser it is possible to locally heat a very small volume of air up to ionization. This process then leads to the creation of a small plasma zone which rapidly expands in all directions. [START_REF] Rossignol | Validation of DLR's sound shielding prediction tool using a novel sound source[END_REF] proposed a sound shielding prediction tool based on a this type of sound source. The lower intrusivity make this technique especially attractive for shielding/installation tests, which typically have to be performed at relativey small scale to reduce costs. In their article, they apply the method to the study of acoustic shielding by a LNA-1 nacelle. More recently, the methodology was applied by [START_REF] Rossignol | Analysis of the noise shielding characteristics of a NACA0012 2D wing[END_REF] to the analysis of the noise shielding characteristics of a NACA0012 2D wing. This methodology is well suited for the study of acoustic shielding. In the case of Green's function estimation, the GFs must be estimated for each of the focus points requiring one session of acquisition for each of them. [START_REF] Bahr | A plasma-based non-intrusive point source for acoustic beamforming applications[END_REF] showed that the process provides an effective in situ method for array calibration both with and without flow and improves the agreement of beamforming results. [START_REF] Fischer | Beamforming in a reverberant environment using numerical and experimental steering vector formulations[END_REF] applied this methodology using a transmitter as a noise source and showed better beamforming results in a hard-walled test-section. However, this solution may not be an option for large configurations and facilities in view of the wind tunnel operating costs and filled order books.

Numerical methods

Most of the numerical methods presented in this subsection concern the study of acoustic shielding. Acoustic monopoles are placed near a geometry and the total diffracted acoustic field is looked for. This enables to estimate, for example, the noise that will be directed towards the ground when a turbofan is installed behind an aircraft wing. Those installation effects can then be used in the aircraft design chain. While those method doesn't concern, strictly speaking, GF estimation, they can be used for their obtention.

To do so, for each of the focus point, a simulation is performed with an acoustic source at that point. The resulting acoustic pressure signal at each microphone is the convolution and the GF can be extracted easily by passing in the Fourier domain.

Asymptotic methods are well adapted at high frequencies where other simulation methods becomes computationally demanding. A method proposed by [START_REF] Candel | Numerical solution of conservation equations arising in linear wave theory: application to aeroacoustics[END_REF] is based on the derivation of a variational system of differential equations for the geodesic elements of the wave front. The elementary cross-section of the wave front is obtained by integration and the principle of conservation of the field invariant directly yields the field amplitude. In addition, suitable jump conditions are derived for treating specular reflexions at solid boundaries. The method is illustrated by specific problems of interest in aeroacoustics. The method solves numerically the ray equations so that an arbitrary

Chapter 1 Methods for the obtention of the acoustic Green's function three-dimensional flow can be treated. The method has difficulties in the determination of the field amplitude in an homogeneous media and in the presence of caustics.

Ray method is also a natural alternative at high frequencies. It is based on Geometric Optics and Diffraction Theory which assumes that all waves are locally plane waves [START_REF] Keller | Geometrical theory of diffraction[END_REF]. There are several advantages of using a geometrical ray-theory approach.

Since the method is geometric, computationally, it is not dependent on the size of the simulation domain nor on the frequency. This method was used by Agarwal et al. (2007) to calculate acoustic shielding by the silent aircraft airframe. The transformation of [START_REF] Taylor | A transformation of the acoustic equation with implications for windtunnel and low-speed flight tests[END_REF] was used to take into account a potential flow around the aircraft. The main drawback of this kind of method is the difficulty to treat complex flow because they assume that rays are straight lines.

The boundary element method BEM or its accelerated variant, the Fast Multipole Method solves the monochromatic version of the wave equation, called the Helmholtz equation [START_REF] Burton | The application of integral equation methods to the numerical solution of some exterior boundary-value problems[END_REF][START_REF] Banerjee | Boundary element methods in engineering science[END_REF][START_REF] Coifman | The fast multipole method for the wave equation: a pedestrian prescription[END_REF]. The problem on the whole domain is reformulated as a problem on the boundaries. The boundaries must be discretized into surface elements and around 10 elements per wavelength are necessary to ensure a good accuracy of the results. This method is well suited for the study of the installation effects at low Mach number [START_REF] Delnevo | Numerical methods: Fast multipole method for shielding effects[END_REF]. It was used by [START_REF] Lummer | Validation of a model for open rotor noise predictions and calculation of shielding effects using a fast bem[END_REF] In this work, "volumic" methods will be used for the simulation. The simulation is carried on the entire volume and a discretization of the whole volume is necessary. While those methods are computationally expensive they can handle a complex geometries and flows. [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationary computational fluid dynamics[END_REF] proposed an approach using excited CFD simulations with system identification (CFD/SI). Thermoacoustic systems are modelled as networks of acoustic multi-ports, where each multi-port corresponds to a certain component of the system, e.g., air or fuel supply, burner, flame, combustor and suitable terminations. In their approach, instationary CFD simulation is performed with broadband excitation at boundary conditions. The system is modeled as a linear time invariant (LTI) system and the application of system identification [START_REF] Ljung | System identification[END_REF] techniques are applied to time series data to extract the unit impulse responses of the multi-port [START_REF] Sovardi | CFD-based modelling of sound generation in ducted discontinuities[END_REF]. [START_REF] Föller | Aeroacoustic characterization of T-junctions based on large eddy simulation and system identification[END_REF] successfully applied the to the aeroacoustic characterization of Tjunctions [START_REF] Föller | Advances in identification techniques for aero-acoustic scattering coefficients from large eddy simulation[END_REF] using large eddy simulation and system identification (LES/SI). We will see in chapter 2 that acoustic GFs can also be seen as the impulse responses of a particular LTI system.

A totally numerical method has been proposed by [START_REF] Kornow | CAA simulation of refraction of sound waves on planar wind tunnel test section with its experimental verification[END_REF] on the use of steady Reynolds averaged NavierStokes (RANS) simulation to describe the time-averaged motion of turbulent flow. Then, in a second step, an acoustic source is placed within the flow and the sound propagation through the flow is computed numerically using Computational AeroAcoustics (CAA). The method was recently applied to the study of trailing-edge noise, slat noise, and jet noise [START_REF] Ewert | CAA broadband noise prediction for aeroacoustic design[END_REF]. However, this study aimed to investigate the effect of the flow on the acoustic propagation and not for the purpose of aeroacoustic imaging.

(p caa m ) m (p exp m ) m (s i ) i (b g i,m ) i,m Numerical#model# CONTENTS 4 l p (x) = ⇧ ⌥ 1 T T ⌦ 0 |x(t)| p dt ⌃ 1/p ( 42
)
s i ⌅ s j (t) = 1 2 NF u=1 a (i,j) u cos ⇥ 2⇥f u t ⇤ (i) u + ⇤ (j) u ⇤ = 1 2 x (i,j) (t) (43) i,n j,n 0 = (s i ⌅ s j ) (t n t n 0 ) ⌅ i, j ⌅ J1, N S K n, n ⌅ J0, N V 1K c i,n = (s i ⌅ p M ) (t n ) ⌅ i ⌅ J1, N S K n ⌅ J0, N V 1K (44) FFT {G M,i } (f ) =? (45)
a u = C te (46) (↵ g i,m ) i = argmin (gi,m) i p caa m i g i,m ⇥s i 2 2 s. t. i ⇧g i,m ⇧ q ⇤ q (47) (↵ g i,m ) m = argmin (gi,m) m p caa i i g i,m ⇥s m 2 2 s. t. m ⇧g i,m ⇧ q ⇤ q (48)
The present PhD dissertation is the continuation of the PhD of [START_REF] Pene | Antennerie numérique pour la caractérisation de sources aéroacoustiques en milieu complexe[END_REF] and lies between the works of [START_REF] Kornow | CAA simulation of refraction of sound waves on planar wind tunnel test section with its experimental verification[END_REF] and [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationary computational fluid dynamics[END_REF]. The procedure of acoustic GF estimation for aeroacoustic beamforming is represented in the figure 1 [START_REF] Ljung | System identification[END_REF] can be applied to extract the acoustic GF ( g i,m ) i,m that are also the impulse response of the system. The estimated GFs can then be used to perform acoustic beamforming using methods described in the section 1.2.2.

It will be shown later in section 2.2.2 that the MISO system identification problem is ill-posed. To regularize it, [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationary computational fluid dynamics[END_REF] added a filter's memory time T M or reverberation time [START_REF] Huang | Acoustic MIMO signal processing[END_REF]. It consists of adding to the model the a priori information that most of the energy of the impulse responses is in this interval of time T M . The number of unknown parameters is reduced and it results in an improvement of the condition number. In the present work, the SI problem will be regularized by adding a sparsity constraint [START_REF] Huang | Acoustic MIMO signal processing[END_REF]. GFs are estimated as a sum of a small number of Dirac delta functions corresponding to acoustic paths from the focus point to the microphone position. The inclusion of the sparsity constraint in the model will also enable to reduce the computational cost.

The mathematical background of the acoustic GF was exposed. It can be formulated as the impulse response of a LTI system representing the medium. The knowledge of the GF makes possible to apply acoustic imaging techniques for complex configurations.

Beamforming based on an approximative GF may lead to errors in characterization of the location and/or amplitude of the sources. Several methods are available for the determination of GFs. Analytical methods consist in the mathematical resolution of the system of PDE but is not applicable on a wide variety of cases. The GF is difficult to obtain especially in the presence of scattering objects and complex flows. Experimental approaches are difficult to apply mainly because it requires extensive human resources, facilities and time. Numerical simulation is an interesting alternative and especially the CAA simulation. It enables to simulate the propagation of sound in the presence of complex geometries and flows but requires a high computational cost. In the following, it is envisaged to use CAA simulation for the estimation of the GFs. A sparsity constraint, valid only for configurations where the GFs can be considered as sparse, is added. This constraint will enable to estimate all the GFs from a single simulation in order to significantly reduce the computational cost.

Chapter 2

Study, regularisation and resolution of the sparse estimation problem

In the previous chapter, a multisource approach was proposed for the estimation the GFs. This approach consists in estimating, for a given acoustic imaging configuration, all the required Green's functions from a single CAA simulation. In this chapter we start by giving the proof that the multisource GFs estimation problem is ill-posed. In order to regularize it, it is necessary to add a priori information on the solution to the model. In the case of external aeroacoustic imaging, in the absence of resonance or reverberation, the acoustic GFs can often be considered sparse. Hence, in that case, the GFs are composed of a small number of spikes. Each spike corresponds to one acoustic path from the focus point to the microphone position. We show in this chapter that the adding of a constraint on the l 0 -or l 1 -norm of the solution enables to regularize the GFs estimation problem. The adaptations of several algorithms for solving this problem are presented and compared.

Simulation of the direct problem

In this section, the methodology for the numerical estimation of the GFs is presented. To obtain the sets of signals needed for the estimation of GFs, the propagation of acoustic waves is modeled by the linearized isentropic Euler equations (1.1). A CAA simulation is performed in order to propagate well-designed signals from potential sources located to a set of focus points through a moving medium. The resulting signals are recorded by virtual microphones. The GFs are finally identified as the transfer functions (or impulse responses) between the source points and the microphones. This step is addressed in the next section.

Computational aeroacoustic simulation

ONERA's code sabrina v0 [START_REF] Redonnet | Numerical simulation of propagation of small perturbations interacting with flows and solid bodies[END_REF] is used to solve the linearized Euler equations in the non-conservative form with disturbances on body fitted structured mesh. The code has been validated on a wide range of complexity up to the simulation of fan interaction noise in 3D coaxial engine [START_REF] Redonnet | Numerical simulation of the downstream fan noise of 3D coaxial engines[END_REF][START_REF] Polacsek | Numerical simulations of fan interaction noise using a hybrid approach[END_REF][START_REF] Redonnet | Computational aeroacoustics of realistic co-axial engines[END_REF][START_REF] Redonnet | Numerical study of acoustic installation effects with a computational aeroacoustics method[END_REF]. This code uses high order finite difference schemes (up to the sixth order), high order filtering schemes (up to the 14th order) and an explicit (third order) time scheme. In the frequency domain, the GFs need to be determined only on a given frequency band [f min , f max ]. Low-dispersive and low-dissipative acoustic propagation is ensured by using at least 10 grid points per wavelength at the highest frequency of interest, f max . To ensure outflow boundary conditions, a combination of a smooth stretching grid and [START_REF] Tam | Radiation and outflow boundary conditions for direct computation of acoustic and flow disturbances in a nonuniform mean flow[END_REF] conditions is used at the boundaries of the domain. To evacuate acoustic waves accurately, the sources of the problem are located at least at a distance of one wavelength from the outflow boundaries (relatively to the smallest frequency of interest, f min ). The propagation of waves outside the frequency range [f min , f max ] may result in higher numerical dissipation and/or spurious reflections at the domain boundaries.

For the i-th source, the forcing term is chosen to be multiplied by a normalized Gaussiandistribution:

c 2 0 ∂Θ i ∂t = 1 3 exp -π x -x i 2 2 × s i (t), (2.1)
where represents the spatial extent of the acoustic source and s i the input excitation signals. The parameter must be taken sufficiently large with respect to the cell size in order to ensure computational stability, while being sufficiently small with respect to the smallest acoustic wavelength to reproduce an acoustically compact noise source. For the simulations performed in the following, = 2∆x where ∆x is the mesh step. The validation of the injection of acoustic monopoles in CAA simulation is presented in the appendix A.

Excitation signals

To

s i (t n ) = N F u=1 a (i) u sin 2πf exc u t n + φ (i) u , (2.2)
where N F is the number of frequencies and

f exc 1 = f min ≤ • • • ≤ f exc N F = f max
are the excited frequencies. In order not to favor any particular source or frequency, the amplitude of all the harmonics is taken unitary a

(i) u = 1 ∀(u, i) ∈ 1, N F × 1, N S .
The signal is thus T -periodic with T = 1/∆f . The sampling period ∆t is constant and the discrete time steps are given by t n = n∆t. The number of time steps for one period of signal is called N T and T = N T ∆t.

The phases φ (i)

u are randomly generated in [0, 2π[ with a uniform distribution in order to ensure decorrelation between all the source signals. Better decorrelation can be ensured by increasing the number of frequencies N F . However, it will also imply a longer period and thus an increase in the computational cost.

Sparsity-based GF estimation problem

Most of the regression algorithms, and especially those studied in the present work, aim to minimize the gradient of the least squares criterion. In this section, the necessity of regularization is highlighted. To do so, the least squares criterion is introduced and the expression of its gradient is derived in the section 2.2.1. In the section 2.2.2, the GF estimation problem is shown to be ill-posed. To regularize it, it envisaged to take into account the sparsity in time of the GFs for external aeroacoustical configurations and several regularization strategies are presented in the section 2.2.3.

Gradient of the least squares criterion

For several sources emitting at the same time, the pressure signal at the microphones is given by (1.6) that represents the propagation of sound in the medium:

p m (t) = N S i=1 (g i,m * s i ) (t) ∀m, t, (2.3)
where * is the convolution product, p m (t) is the pressure signal measured at position x m of microphone m and g i,m (t) will be referred to as the GF between source i and microphone m. It represents the acoustic impulse response of the medium between those two points. For each microphone m independently, GFs between the sources and this microphone appear to be solution of a multi-input, single-output (MISO) system identification problem. Since the minimization process is applied independently for each microphone m, the microphone index m will be omitted for better readability purposes. The Least Squares (LS) criteria is thus defined as:

Q (g) = p - N S i=1 g i * s i 2 2 , (2.4)
where * is the convolution product u * v between two discrete time signals u and v is defined as

(u * v)(t n ) = N T -1 n =0 ū(t n -t n )v(t n ), (2.5)
where ū is the periodic extension of signal u (ū(t nt n ) = u(T + t nt n ) for negative times t n < t n and ū = u elsewhere). Using this notation, the Least Squares (LS) criterion reads as:

Q (g) = p - N S i=1 s i * g i 2 2 = p - N S i=1 A i g i 2 2
.

(2.6)

Here, we introduced matrix A i that reproduce the convolution with source signal s i :

A i n,n = (s i (t n -t n )) n, n ∈ 0, N T -1 , (2.7)
with si is the periodic extension of signal s i .

With these notations, the gradient relatively to GF g i is

∂Q ∂g i = -2A i T p - N S i=1
A i g i .

The matrix A i T is the transpose of matrix A i . It reproduces the cross-correlation with source signal s i :

A i T v = s i ⊗ v with (u ⊗ v) (t n ) = N T -1 n =0 ū(t n -t n )v(t n ), (2.8)
where ⊗ is the cross-correlation product and ū is the periodic extension of signal u. More details on the convolution and cross-correlation products can be found in the book of [START_REF] Gray | Fourier transforms: an introduction for engineers[END_REF] on Fourier transforms. Dropping the multiplicative constant, the expression of the gradient c i,n is

∂Q ∂g i (t n ) ( g) c i,n = +s i ⊗   p - N S j=1 s j * g j   (t n ).
(2.9)

Condition number

Ignoring a multiplying factor, the gradient of the criterion with respect to g i (t n ) writes:

c i,n =   s i ⊗   p - N S j=1 (s j * g j )     (t n ).
(2.10)

A fitting descent algorithm tries to minimize this gradient as it will be seen in the next section. The problem can thus be formulated as finding the functions that cancel the gradient and g is solution of the Wiener filter equation:

  s i ⊗ N S j=1 (s j * g j )   (t n ) = (s i ⊗ p) (t n ) i ∈ 1, N S n ∈ 0, N T -1 . (2.11)
The number of equations here is N S N T , equal to the number of unknowns (every GF at every time step {g i (t n )} i,n ). However, those equations are not linearly independent and the rank r of the system (2.11) is smaller than N S N T . First, all these equations are linear combinations of the equations obtained by simply cancelling the residue:

N S j=1 (s j * g j ) (t n ) = p(t n ) n ∈ 0, N T -1 , (2.12)
meaning that r ≤ N T . Besides, for periodic signals the rank is even lower. The direct discrete Fourier transform F and its inverse F -1 of a signal are respectively given by

U (f k ) = F {u} (f k ) = N T -1 n=0 u(t n )e -i2πtnf k k ∈ 0, N T -1 , (2.13) with f k = k∆f = k/T and u(t n ) = F -1 {U} (t n ) = 1 N T N T -1 k=0 U (f k )e i2πtnf k n ∈ 0, N T -1 . (2.14)
With this expression, and because we are measuring one period of converged signals, the convolution theorem apply. It states that the Fourier transform of the convolution is the pointwise product of the Fourier transforms F {u * v} = F {u} • F {v}. The discrete Fourier transform being a linear operation, applying it to system (2.12) gives, at frequency f k :

N S j=1 S j (f k )G j (f k ) = P (f k ) k ∈ 0, N T -1 , (2.15)
where S j , G j and P are the discrete Fourier transforms of s j , g j and p respectively. The discrete Fourier transform is invertible and, thus, conserve the rank so that system (2.12) has the same rank as system (2.15). As only a small number N F < N T of frequencies are excited during the simulation, most of these equations are always true. The non-trivial equations are those corresponding to the excited frequencies

N S j=1 S j (f exc u )G j (f exc u ) = P (f exc u ) u ∈ 1, N F , (2.16)
Thus we have shown that if multisine signals exciting frequencies {f exc u , 1 ≤ u ≤ N F } are used for the excitation, and assuming that one period of converged signal is recorded, the rank r of system (2.11) is smaller than the rank of system (2.16) that is itself smaller than N F because it contains only N F equations. The problem of GF estimation is thus severely ill-posed (r ≤ N F N S N T ) justifying the need of a regularization term.

An approach for the resolution of (2.

= (s i ⊗ s j ) (t n -t n ) i, j ∈ 1, N S n, n ∈ 0, N T -1 b i (t n ) = (s i ⊗ p) (t n ) i ∈ 1, N S n ∈ 0, N T -1 . (2.18) Auto-correlations (s i ⊗ s i ) (t n -t n
) appear in the diagonal blocks of the matrix. They represent a physical limitation to the estimation of the GFs. Hence, the only excited frequencies can be estimated in the process [START_REF] Schoukens | Survey of excitation signals for FFT based signal analyzers[END_REF]. Considering offdiagonal blocks, it shows that the pairwise signals cross-correlations need to be low.

The difficulty to solve the deconvolution problem also increases with the number of sources N S . The contributions of all the sources are added in Eq. (1.6) making it difficult to extract GFs between the considered microphone and all the focus points.

We observe in practice that the number of frequencies N F required to estimate GFs increases with the number of sources N S .

Sparsity-based regularization

The GF estimation problem is ill-posed. Adding a constraint to the GFs enables to regularize the problem. The resulting inverse problem can then be solved in the time domain if the input and output signals hold enough information. The problem, with l q -regularization, reads as:

g = argmin g p - N S i=1 (s i * g i ) 2 2 s. t. N S i=1 g i q ≤ λ q , (2.19)
where g stands for the concatenation of all the GFs g 1 , • • • , g N S between the microphone under consideration and all the focus points and g is an estimate for g. The parameter λ q can be a predefined constant or can be obtained by using cross-validation [START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF]. The second solution will be used, λ q that yields the smallest cross-validation prediction error is selected. Details on the cross-validation are given in the section 2.4.

The GF estimation problem consists in finding a representation of the pressure at a microphone p m in the dictionary D consisting in the

N S N T signals representing the source signals delayed in time D = {s (n) i } i,n with ∀(i, n) ∈ 1, N S × 0, N T -1 : s (n) i (t n ) = si (t n -t n ) n ∈ 0, N T -1 , (2.20)
where si is the periodic extension of signal s i . With these notations, the GF estimation problem (2.19) can be reformulated as

g = argmin g p - N S i=1 N T -1 n=0 g i (t n )s (n) i 2 2 s. t. i g i q ≤ λ q .
(2.21)

The purpose of linear regression algorithms is to find, from a family of vectors called covariates {u 1 , • • • , u N }, a sparse linear combination that approximates a given vector v.

With α = (α 1 , • • • , α N ) the coefficients of this linear combination, the OMP algorithm solves the problem

α = argmin α v - N i=1 α i u i 2 2 s. t. α q ≤ λ q .
(2.22)

The vector α is the vector of regression coefficients, the family of vectors

{u 1 , • • • , u N }
is the dictionary D and α is an estimate for α.

The similarity between problems (2.21) and (2.22) makes possible the use of linear regression algorithm to solve the GF estimation problem (2.19). However, most of the minimization algorithms require the knowledge of the gradient. Expression of the α LS

α 2 α 1 !"#$%&"'()*+,"& α ridge (a) l 2 -constraint (ridge) α LS α 2 α 1 !"#$%&"'()*+,"& α stepwise (b) l 0 -constraint (stepwise) α LS α 2 α 1 !"#$%&"'()*+,"& α l0.8 (c) l 0.8 -constraint Figure 2
.1: Visualization of a l q -constrained minimization problem with q=2, 0, and 0.8. The iso surfaces of the LS criterion are plotted going from blue (small values) to red (high values) and the constraint region α q ≤1 is plotted in cyan: the region is the abcissa and ordinate axis in the OMP problem.

gradient in (2.22) is expressed by

∂Q ∂α i ( α) u i T   v - N S j=1 α j u j   .
(2.23)

The naive approach consists of computing independently the gradient of the least squares

criterion ∀(i, n) ∈ 1, N S × 0, N T -1 : ∂Q ∂g i (t n ) ( g) s (n) i T   p - N S j=1 N T -1 n =0 g j (t n )s (n ) j   = s (n) i T µ.
(2.24)

The computation of the gradient using this formula has a O N S N 2 T complexity if vector µ is computed initially and then used for the calculation of the gradient. The use of expression (2.10) is preferred in the following. It enables to reduce the algorithmic complexity to O (N S N T log 2 (N T )) by using high speed convolution and correlation [START_REF] Stockham | High-speed convolution and correlation[END_REF]. the regions where α satisfies the constraint (here α 2 ≤ λ 2 ). Because of this penalty, the solutions with small l 2 -norm are privileged. This enables a trade-off between the minimization of the LS criterion Q and the aspect of the solution. This penalization is of prior interest in the case of ill-posed problems because the solution space is an affine subspace and solutions with very high norm are also available.

In general q is taken in [0, 2]. The ridge regression problem corresponds to the case q = 2 α = argmin α v - N i=1 α i u i 2 2 s. t. α 2 ≤ λ 2 , ( 2
A constraint on the l 0 -norm leads to the stepwise regression problem

α = argmin α v - N i=1 α i u i 2 2 s. t. α 0 ≤ λ 0 . (2.26)
It is a variable subset regression problem in the sense that we want to find the solution α that has the smallest number of non-zero components. It has been proved by [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF] to be NP-hard. This property is true for all q < 1 because the constraint region is An increase in the index q of the l q -norm chosen in the minimization problem reduces the sparsity of the solution estimates. Considering the l 1 -norm leads to the lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF])

α = argmin α v - N i=1 α i u i 2 2 s. t. α 1 ≤ λ 1 . (2.27)
This is the smallest value for which the problem is convex making it resolution easier.

The figure 2.2-(a) represents the lasso estimation problem and we can see that the constraint region is a convex polyhedron. An other advantage of constraining the l 1norm is the ability to set coefficients exactly equal to zero because of the sharp corners of the constraint region. This behavior is called "soft-thresholding" and differs from the stepwise regression that operates a "hard-thresholding" as it imposes the number of non zero components with the λ 0 parameter.

α LS α 2 α 1 α lasso !"#$%&"'()*+,"& (a) l 1 -constraint (lasso) α LS α 2 α 1 α l1.2 !"#$%&"'()*+,"& (b) l 1.2 -constraint α LS α 2 α 1 α e-n !"#$%&"'()*+,"& (c) elastic net Figure 2
.2: Visualization of a l q -constrained minimization problem with q=1 or 1.2 and the elastic net problem with t = 0.2. The iso surfaces of the LS criterion are plotted going from blue (small values) to red (high values) and the constraint region α q ≤1

is plotted in cyan.

The lasso can sometimes appear to still be too strict in terms of sparsity of the solution.

The ability of the problem to furnish solutions with numerous coefficients equal to zero is lost when q > 1 (?). For instance, with a constraint on the l 1.2 -norm shown in the figure 2.2-(b), the corners are not sharp and the solution has a few chance to be in the coordinates axes. To remedy this problem, [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] proposed the elastic net problem

α = argmin α v - N i=1 α i u i 2 2 s. t. (1 -t) α 1 + t α 2 ≤ λ, (2.28) with t ∈ [0, 1].
The parameter t enables a trade-off between lasso and ridge regression.

The elastic net problem is represented in the figure 2.2-(c) with t=0.2. Although visually very similar, the elastic-net has sharp corners, while the q=1.2 penalty does not.

In the case of external aeroacoustical configurations, GFs can be considered as sparse.

The GFs are estimated as a sum of time-delayed Dirac delta functions and solutions with the minimum number of spikes are privileged. During the last two decades, there has been an increasing interest on the l 1 -regularization. This has led to numerous algorithms, among them the lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and the elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] algorithms. The l 1 -penalty is popular because it leads to a convex estimation problem that can be solved efficiently using convex optimization methods [START_REF] Friedman | The elements of statistical learning[END_REF]. [START_REF] Efron | Least angle regression[END_REF] proposed the LARS algorithm that can be modified to efficiently solve those regression problems. While the l 0 -penalty hasn't all those properties, the Orthogonal Matching Pursuit algorithm [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] (or forward stepwise regression) enables to produce an approximate path through data. The OMP algorithm is often described as less stable than other l 1 -constrained algorithm that has a more strong mathematical background. At the point that [START_REF] Lin | In defense of l 0[END_REF][START_REF] Lin | A risk ratio comparison of l 0 and l 1 penalized regression[END_REF] have declared themselves defenders of the l 0 -regularization. However OMP algorithm exhibit some interesting properties: a fast convergence and an absence of shrinkage on the estimated coefficients. The practical implementation of the OMP, lasso and elasticnet is described in the section 2.3. Those algorithms will be compared further in the section 2.5 in this chapter.

Adaptation of linear regression algorithms 2.3.1 The Orthogonal Matching Pursuit

The OMP algorithm solves the l 0 -constrained GF estimation problem g = argmin

g m p - N S i=1 (s i * g i ) 2 2 s. t. N S i=1 g i 0 ≤ λ 0 .
(2.29)

It starts with the zero solution [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] and iterates in order to minimize the gradient 2.9. At each iteration, the component that maximizes the gradient is added to the active set (set of all the non-zero components of the GFs) and a linear system is solved in order to determine the amplitudes of these components. The algorithm was originally proposed for the wavelet decomposition in order to compute sparse representations of signals in possibly overcomplete dictionaries. In the case of external aeroacoustics, the algorithm can be used to exploit the sparse characteristics of the GFs. Most of the energy of the GFs is concentrated in a few echoes. It is possible to give a good approximation of the time-domain GFs as a sum of pulses, each pulse corresponding to one acoustic ray between the source and microphone positions.

At each iteration q of the algorithm, one non-zero component is added to one of the GF estimates. The added component is the one which maximizes the gradient (2.10).

A projection on the active set is then performed and all the coefficients are updated in order to decrease the LS criterion (2.4) as much as possible. In this sense, the algorithm behaves like if it is l 0 -constrained. For λ 0 fixed, the different steps of the procedure are described as follows:

Step 0 Initialization of the iteration index q = 0, the algorithm starts with the intercept:

g (q=0) i (t) = 0 ∀i ∀t A (q=0) i = ∅ ∀i A (q=0) = ∅ (2.30)
where g i is the estimate for GF g i . Set A i is the active set relative to source i. It contains all the time steps n where g i (t n ) is not zero. The global active set A contains all the tuples (i, n) for which g i (t n ) is not zero. Superscript (q) is here to differentiate between the value at the previous iteration and the new one.

Step 1 Increment the iteration index q ← q + 1. Compute the gradient using (2.10) based on GFs estimated at the previous step (q -1). Then determine the argument of its maximum in absolute value outside the active set and add it to the belonging active sets:

c (q-1) i,n = +s i ⊗   p - N S j=1 g (q-1) j * s j   (t n ) ∀i ∀n (2.31) ( i (q) , n (q) ) = argmax (i ,n )∈(A (q-1) ) c c (q-1) i ,n
(2.32)

A (q) = A (q-1) ( i (q) , n (q) ) A (q) i =    A (q-1) i n (q) i = i (q) A (q-1) i i = i (q)
(2.33)

Step 2 We can now look for GF estimates as a finite sum of Dirac delta functions:

g (q) i (t) = n∈A (q) i β (q) i,n δ(t -t n ) ∀i (2.34)
The vector of coefficients β (q) = β (q)

i,n (i,n)∈A (q) ∈ R q , representing the amplitudes of the spikes, minimizes the LS criterion:

β (q) = argmin β (q) ∈R q    p(t) - (i,n)∈A (q) β i,n s i (t -t n ) 2    (2.35)
The minimum is obtained for β (q) solution of the system:

Γ A (q) β (q) = α (q) ∈ R q (2.36)
Where: Γ

A (q) = [(s i ⊗ s j ) (t n -t n )] (i,n)∈A (q) (j,n )∈A (q) ∈ R q×q α (q) = [(s i ⊗ p) (t n )] (i,n)∈A (q) ∈ R q
(2.37)

Because we are doing a projection, the solution can completely change in one iteration of the algorithm. It appears that the matrix Γ

A (q)
is a submatrix of the cross-correlation matrix Γ presented in (2.18). The inversion of the linear system at each step can be avoided. Indeed, Γ

A (q-1)
is a Grammian matrix of a set of linearly independent time

signals {s i (t -t n )} (i,n)∈A (q) and Γ A (q)
is a Grammian matrix corresponding to the previous set of vectors increased by one. [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] proposed to use the inverse matrix Γ -1

A (q-1)
at iteration (q -1) for constructing by blocks the inverse matrix at the next step Γ -1 A (q) . It enables to avoid matrix inversion at each step and decrease the algorithmic complexity. Details on this modification are given in the appendix D.

Step 3 Finally GF estimates are updated:

g (q) i (t) = n∈A (q) i β (q) i,n δ(t -t n ) ∀i (2.38) Remaining steps Repeat step 1-3 while i g (q) i 0 = q < λ 0 .
As we will see, the algorithm succeeds well in giving a sparse solution. However, the l 0 -constraint can appear to be too strict in some cases resulting in a bad estimation of the GFs. This is the case when the input and output signals are too short and does not bring enough information for the algorithm to train on and to estimate relevant GFs.

In the worst case, some GFs can even be skipped (estimated as zero functions while it is not the case). It appears in practice that the GFs for focus points located in hidden areas are more likely to be skipped. Considering longer multisine signals with a higher number of harmonics N F will improve the capacity of the algorithm to extract those points.

The LARS-lasso and LARS elastic net algorithms

The elastic net algorithm [START_REF] Efron | Least angle regression[END_REF][START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] is used to solve the GFs estimation problem with constraints on both the l 1 and l 2 -norm of the solution. It deals with correlated estimates but it has nothing to do with correlated source signals.

The method distributes the energy equally between estimates but that may not be exact with GF estimation. The minimization problem reads as:

g m = argmin g m    p m - N S i=1 (g i,m * s i ) 2 2 + λ 2 N S i=1 g i,m 2 2    s. t. i g i,m 1 ≤ λ 1 , (2.39)
with λ 2 a tuning parameter. [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] explain that it is equivalent to solving a problem of the form (2.28) with t = λ 2 /(λ 1 + λ 2 ). The algorithm reduces to the lasso algorithm if λ 2 =0 [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

Like OMP, it is an iterative gradient based algorithm. For λ 1 and λ 2 fixed, the various stage of the procedure are described as follows:

Step 0 Initialization of the iteration index q = 0, the algorithm starts with the intercept:

g (q=0) i (t) = 0 ∀i ∀t A (q=0) i = ∅ ∀i A (q=0) = ∅ (2.40)
where g i is the estimate for GF g i . Set A i is the active set relative to source i. It contains all the time steps n where g i (t n ) is not zero. The global active set A contains all the tuples (i, n) for which g i (t n ) is not zero. Superscript (q) is here to differentiate between the value at the previous iteration and the new one.

Step 1 Increment the iteration index q ← q + 1. Compute the gradient using (2.10) based on GFs estimated at the previous step (q -1). Then determine the argument of its maximum in absolute value outside the active set and add it to the belonging active sets:

c (q-1) i,n = +s i ⊗   p - N S j=1 g (q-1) j * s j   (t n ) -λ 2 g (q-1) i (t n ) ∀i ∀n (2.41) ( i (q) , n (q) ) = argmax (i ,n )∈(A (q-1) ) c c (q-1) i ,n
(2.42)

A (q) = A (q-1) ( i (q) , n (q) ) A (q) i =    A (q-1) i n (q) i = i (q) A (q-1) i i = i (q)
(2.43)

As it will be seen in steps 2 and 3, the absolute value of the gradient on the added component ( i (q) , n (q) ) must be equal to the absolute value of the gradient on the active set A (q-1) . Thus, the absolute value of the gradient on the active set A (q) are all equal to C (q-1) .

Step 2 The direction of equiangular descent is the direction that makes decrease the absolute value of the gradient the same way on all the active set. To find the direction of equiangular descent, consider β (q) the concatenation of the values taken by g (q-1) n on the active set A (q)

β (q) = g (q-1) i,n (τ k ) (i,k)∈A (q) ∈ R N (q) , (2.44)
where N (q) is the cardinal of the active set A (q) . We are looking for a direction of descent ω (q) such that the absolute value of the gradient decreases from C (q-1) to ( C (q-1)γ)

when the Green's function components β (q) become β (q)γω (q) . The equiangular direction of descent is solution of the linear system:

G A (q) ω (q) = (q) ∈ R N (q) , (2.45)
where:

G A (q) = 1 1 + λ 2   (s i ⊗ s j ) (τ k -τ l ) (i,k)∈A (q) (j,l)∈A (q) + λ 2 I N (q)   (q) = sign c (q-1) i,k (i,k)∈A (q)
.

(2.46)

I N (q)
is the identity matrix of size N (q) and sign(x) is the sign of the real scalar x.

The inversion of the linear system at each step can be avoided. Indeed, G

A (q-1)
is a Grammian matrix of a set of linearly independent vectors and G

A (q)
is a Grammian matrix corresponding to the previous set of vectors increased by one. The modification of [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] enables constructing by blocks the matrix G -1

A (q) using G -1 A (q-1)
. Thus, using ω (q) solution of (2.45), ∀γ ≤ C (q-1) it comes:

β (q) -→ β (q) -γω (q)
=⇒ max

(i,k)∈A (q) c (q-1) i,k -→ ( C (q-1) -γ) (2.47)
Moving on this direction enables to decrease the gradient on the active set but may also influence the other components of the gradient. Therefore, γ is chosen as the smallest positive value of γ such that some new index joins the active set (ie. gradient equal to the maximum value).

Step 3 The amplitude of displacement is given by:

γ = min + (i,k)∈(A (q) ) c    C (q-1) -c (q-1) i,k 1 -a (q-1) i,k , C (q-1) + c (q-1) i,k 1 + a (q-1) i,k    (2.48)
Where a

(q-1) i,k
is the cross-correlation of source signals with the direction of descent:

u A (q) (t) = (j,l)∈A (q) ω (q) j,l s j (t -τ l ) a (q-1) i,k = (s i ⊗ u A (q) ) (τ k ) (2.49)
Step 4 Following the Lasso path consists of enforcing the continuous increase of the l 1 norm at each iteration. To do so, non-zero coordinate β (q)

i,k must agree with the sign

(q) i,k of the gradient c (q) i,k : sign β (q) i,k = (q) i,k = sign c (q) i,k
(2.50) [START_REF] Efron | Least angle regression[END_REF] give more details about this modification. To apply the method we need to compute:

γ = min + (i,k)∈A (q)    - β (q) i,k ω (q) i,k    (2.51)
Two situations can then occur: γ can be bigger or smaller than γ. If γ ≥ γ, Step 3bis is finished and no modifications occur on this step. However, if γ ≤ γ, the amplitude of displacement is taken as γ ( γ = γ) and the argument of the minimum in (2.51) is removed from the active set. As one element has been removed from the active set no element is added in the next iteration (Step 1 is skipped). Again, there is a way of constructing by blocks the matrix G -1

A (q) using G -1 A (q-1)
if one vector has been removed from the set of the Grammian matrix.

Step 5 Finally the iteration index is incremented (q = q+1) and the solutions (Green's functions) are updated :

g (q) i,n (τ k ) = g (q-1) i,n (τ k ) + ω (q) i,k ∀(i, k) ∈ A (q) (2.52)
The maximum absolute value of the gradient in (2.41) has been decreased of γ:

C (q) = max ∀i,k ∂R ∂g i,n (τ k ) (g (q) n ) = C (q-1) -γ (2.53)
Remaining steps Repeat step 1-5 while g (q)

n 1 < λ 1 .

Cross-validation stopping criteria

The linear regression algorithms require to set a stopping criterion. Hence, two phases appear during the minimization process: (i) the fitting and (ii) overfitting phases. During the fitting phase, GF estimates are updated in a way that the distance between them and the exact GFs is reduced. The algorithm then reaches an overfitting phase in which the algorithm is trying to extract more information than what the training dataset holds.

In their review on modeling issues, [START_REF] Maier | Neural networks for the prediction and forecasting of water resources variables: a review of modelling issues and applications[END_REF] explain that the overfitting appears due to the high ratio between the number of equations and the number of unknowns. [START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF] proposed the cross-validation to differentiate between the fitting and the overfitting phases. It requires a second set of input signals and pressure at the microphones, the validation set

S = {( s i ) i , ( p m ) m }.
The accuracy of the estimated GFs is calculated as their ability to generalize on the validation dataset that has not been utilized in the training process.

This works as follows. The minimization process (2.19) is applied independently for each microphone index m. We define Q and Q the criteria on the training and validation set:

Q ( g m ) = p m -N S i=1 s i * g i,m 2 2 p m 2 2 and Q ( g m ) = p m -N S i=1 s i * g i,m 2 2 p m 2 2
.

(2.54)

The evolution of this two criteria in function of the iteration number are shown in However the solution differs more outside this interval where numerous spurious spikes have been added to the GFs. The cross-validation appears to give a good estimation of the limit between the fitting and the overfitting phase in this case. The increase For all the cases envisaged in the following, the sample time is taken as ∆t = 10 -6 s.

(i) (ii) (i) (ii) (i) (ii)

Comparison of the convergence of the regression algorithms

In order to compare the various estimation methods, the lasso and elastic net algorithms were also applied to the test case of the previous section. The cross-validation enables to make visible the cap between the fitting and the overfitting phases but suffers from a very slow convergence. For this algorithm, we observe a dropping of the condition number of the matrix Γ

A (q)
before the apparition of the overfitting phase.

The lasso appears to give better GF estimates than the elastic net whatever the value chosen for the λ 2 parameter: the validation curve is always under the one corresponding the elastic net.

The time domain GFs obtained with the lasso at the iteration 11, 200, 300 are presented in the figure 2.6 and differs from those given by the OMP for the same case. After 11 iterations, the OMP algorithm had estimated the main spikes and their amplitude.

The lasso has found the times of arrival of the spikes but the amplitudes given to the coefficients is very low. In fact, the lasso algorithm start by choosing the useful coefficients and then make them grow together. This difference of behavior can also be seen in the figure 2.5 where the slope at the origin for the OMP is very high when compared to the other regression methods. The lasso algorithm is still in the fitting phase after 300 iterations but the convergence is very slow and the GF estimates are (2015) used this regularization for the estimation of the GF in a 2D sparse. In that case, the free-field GF is less sparse and expressed as

g (0,0) (x, t) = - c 0 2π H (t -x /c 0 ) c 2 0 t 2 -x 2 .
(2.55) Thus, the λ 2 parameter represents a lever of action to control the sparsity of the GF.

However, it appears that this spreading of the energy also leads to the apparition of more parasitic spikes.

The comparison of the convergence for a larger number of acoustic sources is shown in the figure 2.8. The same geometry as before is used, the microphone is located at position (0, -2 m, 0) but the number of sources in the segment 0 ≤ x s ≤ 2 m at y=0 m every 0.2 m and z=0 m is now N S =51. The number of frequencies is also increased to N F =701 in order to bring more information to the model. The results are very similar with those of the figure 2.5 only the number of iterations for convergence is higher because more parameters are to be estimated. The minimum of the OMP algorithm is now at iteration 175 instead of 50 for the previous case. The number of focus points N S often reaches several hundreds of as we will see in the chapters 4 and 5. This argument also influenced our choice of the OMP algorithm in the following. In this chapter we proved that the multisource GFs estimation problem is ill-posed and proposed a strategy for its regularization. In the case of external aeroacoustical configurations, GFs can be considered as sparse and this characteristic is used to regularize the problem. The GFs are estimated as a sum of time-delayed Dirac delta functions and solutions with the minimum number of spikes are privileged. The stepwise regression imposes the number of non-zero components of the solution to enforce sparsity and can be solved using the OMP algorithm. The lasso regression solves the problem with a constraint on the l 1 -norm of the solution and can be efficiently solved using convex optimization methods. The sparsity constraint can be relaxed even more with the elastic net algorithm when the GF is less sparse like in the 2D case. We saw that the lasso provides a better estimation of the GFs but has a slow convergence. This is why, in the following, the OMP algorithm will be used for the resolution of the GF estimation problem. The cross-validation stopping criterion and its implementation were also presented.

It prevents overfitting that will lead to a deterioration of the GFs during the estimation process. In the next chapter, the choice of the excitation signals is discussed and in particular alternatives to the choice of randomly-generated phases in the expression (4.2)

of the multisine signals.

Properties of multisine signals

Consider a multisine signal with a defined spectrum density on a given set of frequencies:

x(t) = N F u=1 a u cos (2πf u t + φ u ) , (3.1)
where N F is the number of frequencies and

f 1 ≤ • • • ≤ f N F are the excited frequencies.
All of those frequencies are chosen to be multiple of a frequency f 0 . The signal is thus T -periodic with T = 1/f 0 . The positive real numbers (a u ) 1≤u≤N F define the auto-power spectrum. The root mean square (RMS) value of such a signal is

RMS (x) = 1 T T 0 x(t) 2 dt = 1 2 N F u=1 a 2 u . (3.2)
The RMS value of a signal evaluates the amplitude. The Chebyshev norm (or peak value) of a signal is the maximum of the absolute value of the signal. For a multisine, the Chebyshev norm is defined by its peak value:

x ∞ = max t∈[0,T [ |x(t)| . (3.3)
To finish, the crest-factor (CF) is defined as the ratio between the peak value and the RMS value:

CF (x) = x ∞ RMS (x)
.

(3.4)

The crest-factor of a signal is always larger than 1 and indicates how extreme the peaks are in a waveform. In the following, the source excitation signals are chosen of the form:

s i (t) = N F u=1 a (i) u sin 2πf u t + φ (i) u 1 ≤ i ≤ N S , (3.5)
where

N F , f 1 ≤ • • • ≤ f N F and a (i) u 1≤u≤N F
are auto-power spectrum relative to source i.

In the previous chapters, phases were generated randomly in the interval [0, 2π[. In this chapter, we envisage to construct the set of phases φ

(i) u 1≤u≤N F 1≤i≤N S

Crest-factor minimization

A good overview of the contraints in the choice of excitation can be found in Schoukens

φ u 1≤u≤N F = argmin (φu) 1≤u≤N F x ∞ . (3.6)
As the Chebyshev norm is not differentiable, [START_REF] Guillaume | Crest-factor minimization using nonlinear Chebyshev approximation methods[END_REF] proposed to minimize the continuous l p -norm instead:

l p (x) =   1 T T 0 |x(t)| p dt   1/p . (3.7)
It is then possible to apply gradient based algorithms. The minimization is done for increasing even values of p. The best estimate for the minimization of the l 4 -norm is used to initialize the minimization problem of the l 16 -norm etc. This defines under some regularity conditions an algorithm that converges to the minimax solution (3.6):

φ u u = lim p→∞ argmin (φu) u {l p (x)} . (3.8)
This algorithm enables to reduce the Chebyshev norm of multisine of the form (3.1). This enables to get an excitation signal that maximizes the SNR for given allowable amplitude range. In the following, an adaptation of this algorithm is proposed for the determination of phases that minimizes cross-correlation between signals of the form (3.5).

Cross-correlation minimization

Principles

In this work, we are interested in the minimization of cross-correlation between input signals. Given a set of N S multisine signals s (i) 1≤i≤N S as defined in (3.5), the objective is to find the phases φ

(1)

1 , • • • φ (1) N F , • • • , φ (N S ) 1 , • • • , φ (N S ) N F
that minimize the Chebyshev norm of pairwise signals cross-correlations:

φ (i) u u,i = argmin φ (i) u u,i max i =j s i ⊗ s j ∞ .
(3.9) Given two different input numbers i, j in 1, N S , cross-correlation between the two corresponding signals can be evaluated:

s (i) ⊗ s (j) (t) = 1 T T 0 s i (τ )s j (τ + t) dτ = 1 2 N F u=1 a (i,j) u cos 2πf u t -φ (i) u + φ (j) u = 1 2 x (i,j) (t), (3.10)
where a

(i,j) u = a (j,i) u = a (i) u a (j)
u . The cross-correlation between two such signals appears to be of the form (3.1). Ignoring the multiplying factor it is noted x (i,j) and criteria (3.9) can thus be rewritten:

φ (i) u u,i = argmin φ (i) u u,i max i =j x (i,j) ∞ .
(3.11)

It appears that the minimization problem in this case is close to (3.6) except that minimization has to be done for all the tuples of different inputs. The problem consists of the minimization of multiple signals and the common l p -norm can be considered instead:

l p x (i,j) i =j =   1 T T 0   i =j x (i,j) (t) p   dt   1/p . (3.12)
Again, the minimization is done for increasing values of p. The best estimate for the minimization of the l 4 -norm is used to initialize the minimization problem of the l 16norm etc. This defines under some regularity conditions an algorithm that converges to the minimax solution (3.11):

φ (i) u u,i = lim p→∞ argmin φ (i) u u,i l p x (i,j) i =j . (3.13)
In practice the minimization is not done on all the phases. The phases of the first source can be fixed at their initial values as only difference of phases appear in the cross-correlation expression. This enables to reduce the dimension of the manifold of the solutions and simplifies the problem.

Implementation

In the following, the discrete signal x is the column vector resulting from the concatenation of x (i,j) / 1 ≤ i = j ≤ N S with x (i,j) the column vector of components x

(i,j) n = x (i,j) (t n ) for n ∈ 0, N T -1 . Sampling times are taken as {t n = n∆t, 0 ≤ n ≤ N T -1}
with sampling time ∆t = T /N T . There is then an equivalence between the continuous norm (3.7) and the discrete one:

[l p (x)] p = L p (x) with L p (x) = 1 N T N T -1 n=0 |x n | p . (3.14)
This condition holds if the Shannon sampling condition is respected for signal x p . If not, the discrete norm is only an approximation of the continuous l p -norm.

Minimizing the l p -norm of x consists of minimizing the quantity: (3.15) where e is a column vector of the same size as x whose components are e (i,j) n = x (i,j) n q for i = j, for all time steps n ∈ 0, N T -1 and where q = p/2. Because only even values of p are considered, q is integer and the Jacobian matrix J is then defined as:

L p (x) = 1 N T e T e,
J i,j,n k,u = ∂e (i,j) n ∂φ (k) u =            + q x (i,j) n q-1 a (i,j) u sin 2πf u t n -φ (i) u + φ (j) u if k = i -q x (i,j) n q-1 a (i,j) u sin 2πf u t n -φ (i) u + φ (j) u if k = j 0 else . (3.16)
Using e and J, the solution to the problem can be found using a Gauss-Newton iterative algorithm. The initial phases φ ini are generated randomly in the interval [0, 2π]. They are updated at each step according to:

φ new = φ old -J old T J old + λΛ old -1 J old T e old , (3.17)
with Λ old a positive-definite Levenberg-Marquardt matrix (the diagonal matrix consisting of the diagonal elements of J old T J old ) and λ a parameter that varies with iterations as proposed by [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF]. As the solution improves, λ is decreased and the algorithm approaches the Gauss-Newton method which accelerates convergence to the local minimum. Figure 3.1: (a) Time aspects of signal s 2 with random phases before (---) and after (---) cross-correlation CF minimization. Signal s 1 is unchanged so it is not represented here. (b) Cross-correlation before (---) and after (---) application of the algorithm. The crest-factor is decreased from 2.62 to 1.45 whereas the RMS value is equal to its theoretical value N F /2 = 4.47 and does not depend on the phases.

Cross

Acceleration of the algorithm

The computation of the jacobian matrix can be very slow and needs to be performed at each iteration. Guillaume et al. also give an efficient algorithm, based on fast Fourier transform (FFT), to compute the square matrix J old T J old and the vector J old T e old .

An acceleration using FFT is also possible in our case and details can be found in appendix B.

However, in the present work, J old T J old is a matrix of size N S N F × N S N F . When the number of degrees of freedom (N S N F ) increases, the construction of this matrix becomes very expensive both in terms of computation and memory storage. Moreover, the scalability of the matrix inversion with N S N F is even worse. It appears in practice that it is possible to consider only its diagonal Λ old multiplied by a factor α (α < 1):

φ new = φ old -α Λ old -1 J old T e old .
(3.18)

The algorithm then reduces to a gradient descent (GD) method. The memory storage required for the computation of J T J drops from (N S N F ) 2 to N S N F and its inversion becomes trivial. The parameter α has to be chosen small enough to ensure stability.

In practice taking α = 0.25 was enough for the cases considered in the following. For cases with larger gradients, α may have to be taken smaller. For two source numbers i = j, gradients of s (i) ⊗ s (j) p (t) with respect to phases increase with the index p of the considered l p -norm and with the maximum frequency of s (i) ⊗ s (j) requiring smaller value of α. Ten sources are considered with flat auto-power spectrum containing 100 frequencies f u = 1, • • • , 100 Hz. The algorithm is stopped when the crest-factor variation between two iterations becomes smaller than = 10 -4 two times in a row. These curves were computed on the same machine. Vertical dotted lines stand for 10 iterations of the algorithm and vertical dashed lines delimit the index p of the minimized norm l p .

The number of iterations is increased especially at the end of the convergence. This is not a big counterpart as the main part of the minimization is done in the first iterations.

In Fig. 3.2, the algorithm was applied to the case of N S = 10 inputs with flat auto-power spectrum containing 100 excited frequencies f u = 1, • • • , 100 Hz (a

(i) u = 1 ∀i, u)
. The initial phases are randomly generated in [0, 2π[ giving a crest-factor equal to 3.9692. The Levenberg-Marquardt (LM) algorithm reduces it slightly more (1.83 in 210 iterations) than the GD algorithm (1.87 in 549 iterations). On the other hand, in terms of execution time, the GD algorithm converges in only 42 s. This is way smaller than the 1045 s required by the LM algorithm.

The limitation of these algorithms in terms of computation and memory can be identified. For the LM algorithm, the limitation is the computation and the storage of the J T J matrix of size N S N F × N S N F . Besides the performance of the GD algorithm is limited due to the computation and storage of the cross-correlation vector x of length N T N S (N S -1)/2. In practice, (N S , N F ) = (30, 300) was the limit for the LM algorithm and, using the GD algorithm, it was possible to get to (N S , N F ) = (150, 1500).

Eventually, for the GD algorithm, correlations can be computed on the fly. This enables to reduce the memory storage at the cost of a loss of modularity.

Influence of the number of signals and frequencies

In the previous case only two signals were considered. In certain situations, a large number of input signals may be required and their pairwise cross-correlations have to show a minimum crest-factor. Hence, the RMS value of the cross-correlation is equal to N F /2 whereas the maximum auto-correlation (diagonal elements of the cross-correlation matrix in (2.17) is equal to N F /2. Thus, the ratio between the largest non-diagonal element and the diagonal ones is CF × 2/N F and we want it to be as small as possible. One way of minimizing this ratio is to increase the number of frequencies. This means adding more information in the model as more samples N T may be required. The other way is to find the phases solution of problem (3.11) that minimize the maximum cross-correlation crest-factor. This is what we are interested in in this part.

In The study of the numbers of degrees of freedom and constraints enables to explain this limitation. Hence, the number of degrees of freedom is N dof = (N S -1)N F corresponding to the number of phases to determine. On the other hand, the number of constraints is N cons = N T N S (N S -1)/2, the number of time steps where signals have to be minimized.

The minimization is better when the ratio r = N dof /N cons is large. As r = 1/(2N S )

here, this ratio tends to drop when the number of signals increases making it difficult to minimize the cross-correlation as can be observed in Fig. 3.3. For larger values of N S , Chebyshev norm of the cross-correlation cannot be minimized as much as in the case of

N S = 2.
A cross-correlation crest-factor equal to 1 appears to be the infimum. However, if an a priori value of the delay T M under which the minimization has to be done is given, we can expect to get below this limit. Calling N M the number of time delays for which |t n | ≤ T M , the number of constraints is then expressed as N cons = N M N S (N S -1)/2.

Then N M can now be smaller than the number of frequencies N F . If N M is chosen small enough, a good minimization of the cross-correlation can be obtained even for big numbers of signals. The next section deals with cross-correlation minimization inside a limited time interval.

Adaptation for low-correlation zone signals

For FIR systems the impulse response has a finite duration time also called filter's "memory". Thus, its coefficients h

(i) n = h (i) (t n
) can be neglected for indices n larger than a value N M where N M ∆t = T M has to be larger than the filter's "memory" of previous signal states [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationary computational fluid dynamics[END_REF]). If an upper bound of this caracteristic time is known a priori, then h(t) = 0 for t > T M . Let's denote h the concatenation of the impulse responses h = h (i) n for all the sources i ∈ 1, N S in the considered time range n ∈ 0, N M . The discrete version of Wiener filter equation is still of the form (2.17).

However, the correlation matrix Γ and vector c are now expressed as:

Γ i,n j,n = s (i) ⊗ s (j) (t n -t n ) i, j ∈ 1, N S n, n ∈ 0, N M c i,n = s (i) ⊗ y (t n ) i ∈ 1, N S n ∈ 0, N M . (3.19)
In that case, the problem is driven by the cross-correlation between input signals for delays smaller or equal to T M . By applying the minimization algorithm in this time range only, the Chebyshev norm can be decreased more and can even be smaller than the RMS value taken on all the period. In the algorithm, this is done by taking the Signal s 1 is unchanged so it is not represented here. (b) Crosscorrelation before (---) and after (---) application of the algorithm. The chosen region for low-correlation is represented in gray, here T M = 0.2 s. Cross-correlation crest-factor in the minimization zone is decreased from 2.62 to 0.12.

components of vector e as e (i,j)

n =    x (i,j) n q if t n ∈ [0, T M ] [T -T M , T [ 0 else . ( 3
(i) ⊗ s (j) (-t) = s (i) ⊗ s (j) (T -t)
as periodic signals are considered.

The algorithm can be used for the generation of low correlation zone signals using (3.20).

In Fig. 3.4, the algorithm was applied to the same case as in section 3.3.2 but a minimization zone is considered. This time, the minimization of the cross-correlation is done for small delays |t| ≤ T M where T M = 0.2 s. Periodic input signals verify

s (i) ⊗ s (j) (-t) = s (i) ⊗ s (j) (T -t)
meaning that minimization also has to be done on time interval [T -T M , T [. Again, the time aspect of signal s 2 is slightly modified. The crest-factor of the cross-correlation in the minimization zone is decreased from 2.62 to 0.12. The RMS value of the cross-correlation does not depend on the phases and stays unchanged during the minimization process. However, in this case the energy of the cross-correlation signal tends to move out of the low correlation zone. This enables to get a local crest-factor lower than 1. The local crest-factor is minimized more for small values of T M .

Influence on Green's function estimation

The proposed algorithm is intended to be used for the generation of the excitation signals in the CAA simulation. We remind the expression of the cross-correlation matrix in (2.17) of the estimation problem:

Γ i,n j,n = (s i ⊗ s j ) (t n -t n ) ∀i, j ∈ 1, N S ∀n, n ∈ 0, N T -1 .
This matrix can be decomposed into blocks for each pair of source numbers (i, j)

Γ =          Γ 11 Γ 12 • • • Γ 1N S Γ 21 Γ 22 • • • Γ 2N S . . . . . . . . . . . . Γ N S 1 Γ N S 2 • • • Γ N S N S          , (3.21) with Γ ij n,n = (s i ⊗ s j ) (t n -t n ) ∀n, n ∈ 0, N T -1 . (3.22)
The diagonal blocks do not depend on the phases. Hence, phases disappear in the calculation of cross-correlation (3.10) when i = j. This represents a physical limitation in the estimation process because the only excited frequencies are estimated in the process. The use of multisine signals with minimized cross-correlation enables to reduce the Chebyshev norm of the elements of the non-diagonal blocks or, in other words, the maximum value of the coefficients. The minimization is interesting but in practice the influence on the estimation process is not significant. Estimations were performed for several frequencies and results may even be worse with a set of phases resulting from cross-correlation minimization.

|| • || 1 RMS || • || ∞ CF s sqr A A A 1 s sin 2A √ 2 π ≈ 0.90A A A √ 2 √ 2 s tri A √ 3 2 ≈ 0.87A A A √ 3 √ 3
Table 3.1: Norms of the square, sinus and triangle signals.

Several arguments may explain the poor results obtained. In this chapter, we made the choice of the minimization of the Chebyshev norm. However, it can be shown analytically that the RMS value of the cross-correlation also doesn't depend on the phases of the

s sqr =(-1) k A and k = 2t T s sin =A √ 2 sin (2πt/T ) s tri =(-1) k 4A √ 3 T t - kT 2 and k = 2t T + 1 2 . (3.23)
The value of the l 1 -norm, RMS, Tchebyshev norm and crest-factor of those three signals are listed in table 3.1. It shows that, with a fixed RMS value, signals with the smallest crest-factor have a greater l 1 -norm.

We have seen in the previous section that, if the filter's memory time T M is known a priori, only the part corresponding to small delays is used during the GF estimation process. When this memory time is known and is very small compared to the period T , the proposed minimization method is of sure interest and enables to reduce significantly the cross-correlation in the zone of interest as seen in the figure 3.4. However without this information, method like the pseudo inverse [START_REF] Golub | Calculating the singular values and pseudo-inverse of a matrix[END_REF] or singular value decomposition [START_REF] Golub | Singular value decomposition and least squares solutions[END_REF] may use the whole cross-correlation matrix.

On the other hand, sparse estimation algorithms presented in the previous chapters use the restriction of this matrix on the active set (other elements may however be used for the determination of this active set...). It is then difficult to justify the choice of prioritizing the minimization of the Chebyshev norm over the l 1 -norm.

A minimization algorithm for the generation of multisine signals with low pairwise crosscorrelation has been proposed. The algorithm allows the user to define the auto-power spectrum of the sources and iterates on the phases in order to minimize the common l p -norm with increasing values of p. The minimization is harder for large numbers of input signals. However, in every studied case, the crest-factor of the cross-correlation has been decreased by at least a factor 2. For low memory length systems, an adaptation of the algorithm to the generation of low-correlation zone signals has been presented.

The results show that introducing an a priori information of the time interval on which cross-correlation has to be minimized enables to decrease it more effectively. We saw that the minimization of the Chebyshev norm also results in the increase of the l 1 -norm making questionable the benefits brought by the method. In particular, in the case of the sparse linear regression of interest in the present work, the contributions of this kind of signals are not significant. Consequently, multisine (3.1) with random phases are used for the excitation signals in the following chapters. In the next chapter, the GF estimation procedure is applied to numerical test cases about acoustic imaging in the presenc of a sphere or a NACA0012 2D wing. Estimated GFs will be compared to ones obtained experimentally but no experimental microphone array data will be used until chapter 5.

Chapter 4

Application of the method on numerical test cases

In previous chapters, a methodology has been proposed for the estimation of the Green's function. The methodology enables to estimate all the Green's functions required to perform acoustic imaging. However, we have seen that it is restricted to the framework of extern aeroacoustics because it takes advantage of the sparsity of the Green's functions.

In this chapter, two test cases representative for extern shapes for fuselage and wings are considered. The first case that we address in section 4.1 is the diffraction by a rigid sphere. The acoustic source is located at the opposite side of a microphone line array. The pressure field resulting from the diffraction of a monopole source by a rigid sphere can be derived analytically. In order to prove the advantages of CAA in terms of GF estimation, we will compare CAA estimated GFs with both experimental and analytical ones. The second test case is the diffraction by a NACA0012 wing profile presented in section 4.2. For this test case, the solution is not known analytically and the GFs obtained using the proposed method will be compared with GFs obtained experimentally. Beamforming results will also be presented and we will see that the method can be applied to real industrial configurations with several thousands of focus points. Results presented in this chapter were published in the AIAA Journal [START_REF] Bousabaa | Sparse green's functions estimation using orthogonal matching pursuit: Application to aeroacoustic beamforming[END_REF].

Numerical setup

The simulations were performed with ONERA's sAbrinA-v0 code on a body-fitted mesh. The useful part of the mesh is the zone such that |x| ≤ 0.26 m, |y| ≤ 0.16 m and -0.36 m ≤ z ≤ 0.38 m. The origin of the spatial system is taken as the center of the sphere of radius a = 0.06 m. In the region of uniform mesh, the spatial step is 2 mm. It is small enough to ensure more than 10 points per wavelentgh for the maximum frequency f max = 11 kHz (leading to the minimum wavelength λ min 3.1 cm). A 3% stretching is added at the limits of the domain and combined with [START_REF] Tam | Radiation and outflow boundary conditions for direct computation of acoustic and flow disturbances in a nonuniform mean flow[END_REF] outflow boundary conditions (BC) to evacuate outgoing waves properly. The simulation is led in a medium at rest and in the presence of a potential flow of v∞ = v∞ u x with v ∞ = 45 m/s the flow far from the sphere and u x the unit vector in the x-axis direction (Mach number M = 0.13). This corresponds to the configuration of the experiments described in the previous subsection. ! " #$%&'(#

!"#$ !"#$ !"%%$ !"%$ !"&%$ !"&$ !"'%$ !"'$ !"(%$ !"($ !")%$ !")$ !"!%$
!"# $$# $"# %$# %"# &$# &"# '$# '"# ($# ("# $# !$# (a) x-axis velocity vx ! " #$%&'(# !"# !$# %"# %$# "# $# &"# &%$# &%"# &!$# &!"# &'$# '$# (b) z-axis velocity vz !"#$ !"#$ !"%%$ !"%$ !"&%$ !"&$ !"'%$ !"'$ !"(%$ !"($ !")%$ !")$ !"!%$

Estimation of the Green's functions

In order to validate the CAA simulation approach, we start exciting only one monopole source. This enables to evaluate possible numerical errors and to quantify the quality of the GFs that can be extracted from the simulation. The monopole source is placed at In the following, the focus is put on 3 microphone positions:

• In the line of sight x m = -0.2 m: the source is visible by the microphone. Two main rays are to be considered: one direct ray and a ray resulting from the reflexion on the surface of the sphere. The delay between those two rays will imply positive or negative interferences on the considered frequency range.

• In the shadow zone x m = -0.07 m: the source is hidden by the sphere. The main part of the energy is held by two main rays creeping on the surface clockwisely and counter-clockwisely. In the frequency range [4 kHz,11 kHz], interference between those two rays is mainly destructive and the masking effect is maximum.

• Opposite to the source x m = 0 m: the source, the microphone and the center of the sphere are aligned. Creeping rays come from all the perimeter of the sphere and interfere in a constructive way creating an important increase of the amplitude. The pressure field resulting from the diffraction of a monopole by a rigid sphere can be expressed analytically by [START_REF] Bowman | Electromagnetic and acoustic scattering by simple shapes[END_REF] For the three microphone positions, a good agreement is obtained between analytical, numerical and experimental GFs especially on the phase. In the absence of flow, the GFs extracted from the simulation should equal the analytical solutions. However, because of a small mesh dissipation effect at high frequencies, we notice a small offset in amplitude that reaches 0.5 dB at 11 kHz. This error is negliged here but could be reduced more if a finer mesh is used. However, the consequence will be an increase of the computational cost. The error on the modulus for the experimental GFs is smaller than 2 dB on almost all the chosen frequency range.

!"#$%&'# !" #$%" #$&" #$'" #$(" #" )#$(" )#$'" )#$&" #$%" )!" (#$)'# *#$)'# #$*" #$(" #$!" #" )#$!" )#$(" )#$*" )#$(" )#$!" #" #$!" #$(" !"#$%&'# !" #$%" #$&" #$'" #$(" #" )#$(" )#$'" )#$&" #$%" )!" (#$)'# *#$)'# #$*" #$(" #$!" #" )#$!" )#$(" )#$*" )#$(" )#$!" #" #$!" #$("
The same procedure is followed in the presence of a 45 m/s mean flow (Mach number M = 0.13). Again, several approximations are to be taken into account. In the simulation and for the analytical solution, a potential mean flow is considered. This may induce a difference with the experimental flow especially because of the presence . This shielding effect is underestimated when using the analytical solutions. Finally, for position x m = 0 m, the flow induces a small increase of the amplitude visible on both the experimental and numerical GFs while the flow correction states that the flow has no influence onto the GF for this position.

Thus, the analytical formulation gives an accurate estimation of the GF in the absence of flow. However, as it is shown in this section, the main difficulty is on reproducing correctly the impact of a complex mean flow. The CAA approach enables to take into account complex flows. However, because it is a volumic method, this generally results in a high computational cost that limits its intensive use in a study. In the next section, we will show how the estimation of all the GFs from a single CAA simulation enables to reduce computational cost.

Acoustic imaging

In this subsection, we propose to use the multisource approach for the estimation of all the GFs from one main simulation by solving the GFs estimation problem presented in chapter 2. The whole procedure followed in this section is described in the figure 4.9.

Acoustic GFs are estimated using a CAA simulation and are used to solve a problem of imaging with the classical beamforming method. Because no microphone array data is available, the aeroacoustic measurements are replaced by a CAA simulation to extract a synthetic cross-correlation matrix. The estimated GFs can finally be used to perform acoustic beamforming.

The setup is the same as in the previous section but only the number of sources and microphones and their positions are changed. The 41 microphones representing the sensor array are uniformly placed every 10 mm at z m = -0.3 m and -0.2 m < x m < 0.2 m.

1.#Aeroacous,c#measurements# Physical#configura,on# Numerical#model# Post#processing# microphone# signals############ cross:correla,on# matrix# 2.#Numerical#simula,on# Computa,onal#AeroAcous,c# simula,on# 3.#Green's#func,on#es,ma,on# resolu,on#of#the#GF#es,ma,on#problem#########:# es,mated#Green's# func,ons# (p exp m ) m (b g i,m ) i,m Numerical#model# 4.#Acous,c#beamforming# acous,c#map# synthe,c#source# excita,on#signal# CAA#simula,on# excita,on#signals#at# focus#points# (s i ) i pressure#signals#at# microphones#(p caa m ) m 8m CONTENTS 4 l p (x) = ⇧ ⌥ 1 T T ⌦ 0 |x(t)| p dt ⌃ 1/p (42) s i ⌅ s j (t) = 1 2 NF u=1 a (i,j) u cos ⇥ 2⇥f u t ⇤ (i) u + ⇤ (j) u ⇤ = 1 2 x (i,j) (t) (43) i,n j,n 0 = (s i ⌅ s j ) (t n t n 0 ) ⌅ i, j ⌅ J1, N S K n, n ⌅ J0, N V 1K c i,n = (s i ⌅ p M ) (t n ) ⌅ i ⌅ J1, N S K n ⌅ J0, N V 1K (44) FFT {G M,i } (f ) =? (45)
a u = C te (46) (↵ g i,m ) i = argmin (gi,m) i p caa m i g i,m ⇥s i 2 2 s. t. i ⇧g i,m ⇧ q ⇤ q (47) (↵ g i,m ) m = argmin (gi,m) m p caa i m g i,m ⇥s m 2 2 s. t. m ⇧g i,m ⇧ q ⇤ q (48) (↵ g i,m ) i = argmin (gi,m) i p caa m i g i,m ⇥s i 2 2 s. t. i ⇧g i,m ⇧ 0 ⇤ 0 (49) (↵ g i,m ) m = argmin (gi,m) m p caa i m g i,m ⇥s m 2 2 s. t. m ⇧g i,m ⇧ 0 ⇤ 0 (50)
Figure 4.9: Procedure for acoustic beamforming based on CAA estimated GFs for the application on the sphere test case. During the simulation, focus points are excited using acoustic monopoles with excitation signals (s i ) i and the pressure fluctuation signals (p caa m ) m at the microphones positions are recorded. The GFs are estimated by applying LTI system identification and adding a constraint on the l 0 -norm of the solution. Because no microphone array data is available, the aeroacoustic measurements are replaced by a CAA simulation to extract a synthetic cross-correlation matrix. The estimated GFs can finally be used to perform acoustic beamforming.

On the other side of the sphere, the focalization area is divided into 101 focus points uniformly placed every 4 mm at z m = +0.32 m and for x s in the range [-0.2 m, 0.

s i (t) = N F u=1 cos (2πf u t + φ u,i ) 1 ≤ i ≤ N S , (4.2)
where f u = f min + (u -1)∆f and the phases φ u,i are randomly generated in [0, 2π[ in order to ensure decorrelation between all the source signals. After a transient regime of t trans 4 ms, one period T = 200 ms of signal is recorded. For the validation set S, the sphere. However, we see that the beamforming patterns obtained using the analytical and estimated GFs are of very good quality. Those GFs enable to identify a source masked by the sphere both in terms of position and amplitude. The correct level of 0 dB is found at the position of the source. We can also notice that the beamforming patterns are affected by the presence of the sphere. Its presence increases the secondary lobe amplitudes but reduces the width of the main lobe. At 10 kHz, the maximum side-lobe amplitude are increased to -6 dB when the source is masked by the sphere whereas it only reaches -13 dB when the source is not masked. On the other hand, the width of the main lobe below the sphere is 0.04 m, narrower than in the absence of the sphere (0.064 m).

!"#$%&'# !"# $"# %"# &"# '"# "# ('"# (&"# (%"# ($"# (!"# (#$)'# *#$)'# ")%# ")&# ")'# "# (")'# (")&# (")%# (")&# (")'# "# ")'# ")&#
Finally, the application of the proposed method enabled to reduce the computational cost. Hence, supposing that the GFs are required every 125 Hz, the estimation of the GF for one focus point will require t trans + T = 12 ms. Because one simulation has to be done for each of the N S = 101 focus points, the total simulation time required to estimate the GFs using the monosource approach is t mono = N S (t trans + T ) = 1.12 s. The proposed method needs only t multi = 0.216 s of CAA simulation i.e. 18% of the initial computational cost. If the GFs samples are required every 50 Hz, only 9% of the 2.42 s original computational time is required. Concerning the resolution of the estimation problem (2.29), the cost is negligible in comparison with the simulation. The CAA simulation required around 600 hours on 480 CPU whereas the deconvolution process required 24 hours of Matlab computation without parallel implementation.

Beamforming in presence of a diffracting NACA0012 profile

We now adress the case of the 2D acoustic imaging of a NACA0012 wing profile. For this test case, the solution is not known analytically and the GFs obtained using the proposed method will be compared with GFs obtained experimentally. The spatial sampling of the surfaces results in a high number of focus points. However, the use of the reverse-flow reciprocity will enable applying the method.

Experimental setup

Again, the data were recorded in the DLR AWB anechoic wind tunnel and the experimental setup is presented in figure 4.12. The center of the 0.2 m chord NACA0012 profile wing is taken as the origin of the spatial system. The wing is held from its extremities and is supposed infinite in the Y direction. Again, ONERA's SPARC source, is used for the excitation and the same procedure as the one described in 4.1.1 is used to improve the signal-to-noise ratio. The measurements are performed at a flow speed of 55 m/s.

!"#$ !"#$ !"%%$ !"%$ !"&%$ !"&$ !"'%$ !"'$ !"(%$ !"($ !")%$ !")$ !"!%$ !"#%$

Numerical setup

The simulations were performed with ONERA's sAbrinA-v0 code on a body-fitted mesh. The useful part of the mesh is the zone such that -0.26 m ≤ x, y ≤ 0.26 m and -0.16 m ≤ z ≤ 0.46 m. The origin of the spatial system is at the middle of the 0.2 m chord of the NACA0012 2D wing. In the region of uniform mesh, the spatial step is 2 mm. It is small enough to ensure more than 10 points per wavelentgh for the maximum frequency f max = 11 kHz (leading to the minimum wavelength λ min 3.1 cm).

A 3% stretching is added at the limits of the domain and combined with [START_REF] Tam | Radiation and outflow boundary conditions for direct computation of acoustic and flow disturbances in a nonuniform mean flow[END_REF] outflow BC to evacuate outgoing waves properly. The global mesh contains 78.4 10 6 points and the CFL on the whole mesh is plotted in figure 4.13.

The simulation is led in the presence of a potential flow of 55 m/s in order to be compared with the experiments described in the previous section. This time, the potential mean

Estimation of the Green's functions

In order to have a reference solution for the CAA simulation approach, we start exciting only one monopole source at various focus points. This enables to evaluate possible numerical errors and to quantify the quality of GFs that can be extracted from the multisource simulation. The monopole source is excited using a multisine signal exciting the frequency range [4 kHz, 11 kHz] every ∆f = 125 Hz. The period of the signal extracted from the simulation is T = 8 ms after a transitory time of around 3 ms.

Phases were generated randomly in [0, 2π[. The pressures at the microphone positions are recorded and the GFs are extracted dividing the microphone spectrum by the source

! " #$%&'(# !"# $$# $"# %$# %"# &$# &"# '$# '"# ($# ("# $#
(a) x-axis velocity vx spectrum at the excited frequencies

! " #$%&'(# !"# !$# "# $# %"# %!$# %!"# (b) z-axis velocity vz !"#$ !"#$ !"%%$ !"%$ !"&%$ !"&$ !"'%$ !"'$ !"(%$ !"($ !")%$ !")$ !"!%$ !"#%$
G i,m (f exc u ) = P m (f exc u ) S i (f exc u ) ∀u ∈ 1, N F . (4.3)
This expression holds when only one source i is excited during the simulation and this is the case in this subsection. When the number of sources excited is bigger, the equation (4.3) doesn't hold and the sound propagation model is described by (2.16).

Comparisons between the experimental GFs and the ones obtained using CAA simula- The GFs are obtained with a good accuracy even when the source is hidden below the profile. A shift in phase appears in every case that can reach π at high frequencies.

tion
X Z !"#$ !%#$ !&#$ !"#$ !%#$ !&#$ (a)
This can be due to errors in the post processing of experimental data or to differences between the experimental configuration and the numerical model. The more important intrusivity of the SPARC source in the flow which is placed only at 0.1 m from the diffracting object may also have an influence.

Acoustic imaging

In this subsection, we consider doing aeroacoustic imaging in the presence of a NACA0012

profile wing where all the GFs are estimated from one main simulation. Again, it is a (a) It is included in the plane z = 0.4 m and is at a distance of 0.4 m from the wing. Microphones are uniformly distributed every 0.01 m and the two braces of the array are parallel to x and y-axis.

!"#$%&'# !""# $%"# $!"# &"# '"# "# ('"# (&"# ($!"# ($%"# (!""# (b)

Reverse-flow reciprocity

This time, the number of focus points is very large and solving the estimation problem (2.29) for N S = 4653 is not the best strategy. Hence, in this case the estimation problem is more difficult to solve because it requires a higher number of frequencies due to the fact that more parameters must be estimated for one minimization. The reciprocity principle states that the acoustic response of a medium does not change when the source and microphone are swapped [START_REF] Heaslet | Reciprocity relations in aerodynamics[END_REF]). This principle doesn't hold directly to the case of a moving medium encoutered in aeroacoustics. However, it was pointed out by [START_REF] Lamb | On reciprocal theorems in dynamics[END_REF] that inverting the direction of the steady flow may restore reciprocity [START_REF] Deneuve | A coupled timereversal/complex differentiation method for aeroacoustic sensitivity analysis: towards a source detection procedure[END_REF][START_REF] Padois | Experimental localization of an acoustic sound source in a wind-tunnel flow by using a numerical time-reversal technique[END_REF].

In this section, we propose to assess the applicability of this principle in the case of a

NACA0012 wing placed in a potential flow. A numeric validation of this principle is shown in figure 4.17 in the presence of a NACA0012 diffracting wing. In this example, the estimation of the GF between a source point located near the trailing edge and a microphone is presented. The direct problem consists of exciting the source point using a multisine signal and to measure one period of the resulting pressure signal at the position of the microphone (Figure 4.17-(a)). The frequency domain GF between the two points is obtained by computing the ratio of the Fourier transforms of those signals.

In the reverse approach, shown in figure 4.17-(b), the flow is inverted and the roles of the microphone and the source point are exchanged.

The time aspect of the recorded signals is shown in the figure 4.18 where one period of converged signals is recorded. The aspects of signals obtained with the direct and reverse approach are similar. There appears to be a higher dissipation for the direct approach for which the peaks have a slightly smaller amplitude. The GF in the frequency domain obtained with the direct and reversed approach are shown in figure 4.17-(c) and compared with experimental ones. The direct and reversed approaches give the same GF and are in good agreement with the experiments.

In the present work, this principles enables to inject monopole sources from the microphones to the focus points. The problem of Eq. (2.29) is now solved for each focus point independently and reads as: where N M is the number of microphones and g i now stands for the concatenation of all the GFs g i,1 , • • • , g i,N M between the focus point i and all the microphones.

g i = argmin g i p i - N M m=1 (g i,m * s m ) 2 2 s. t. i g i,m 0 ≤ λ 0 , (4.4) !"#$%" &"#$%" !"#$ !"%$ !"&$ !"'$ !$ (!"'$ (!"'$ !"&$ !$ !"'$ (!"&$ (a) !"#$%" &"#$%" !"#$ !"%$ !"&$ !"'$ !$ (!"'$ (!"'$ !"&$ !$ !"'$ (!"&$ (b)
The procedure followed in this section is described in the figure 4 The GFs are estimated by applying LTI system identification and adding a constraint on the l 0 -norm of the solution. Because no microphone array data is available, the aeroacoustic measurements are replaced by a CAA simulation to extract a synthetic cross-correlation matrix. The estimated GFs can finally be used to perform acoustic beamforming.

u 0 := u 0 (s m ) m 8i (p caa i ) i synthe,c#source# excita,on#signal# CAA#simula,on# CONTENTS 4 lp (x) = ⇧ ⌥ 1 T T ⌦ 0 |x(t)| p dt ⌃ 1/p (42) si ⌅ sj(t) = 1 2 NF u=1 a (i,j) u cos ⇥ 2⇥fut ⇤ (i) u + ⇤ (j) u ⇤ = 1 2 x (i,j) (t) (43) i,n j,n 0 = (si ⌅ sj) (tn tn0 ) ⌅ i, j ⌅ J1, NSK n, n ⌅ J0, NV 1K ci,n = (si ⌅ pM ) (tn) ⌅ i ⌅ J1, NSK n ⌅ J0, NV 1K (44
The reciprocity principle is of particular interest in the case of real aeroacoustic imaging applications for two reasons. On one hand, the number of focus points can reach several thousands whereas the number of microphones in the array is usally smaller or around one hundred. The use of the reciprocity conducts to the problem of the estimation of N M GFs simultaneously instead of N S for the direct problem. Less parameters are to be estimated which results in a reduction of the number of frequencies N F required, which also drives the quantity of information given to the model. This also implies a reduction of the period T of the signal and thus a reduction of the CAA simulation time. On the other hand, reciprocity also enables to extract GFs between microphones and focus points that may be located on the surface of diffracting objects because only pressure perturbations need to be recorded at the focus points. This enables to localize the acoustic sources resulting from the interaction of turbulent flows with structures.

The reverse-flow reciprocity was applied here and a 3D view of the performed CAA which is not true with the flow-corrected free-field GFs. For focus points outside the surface of the wing (x < -0.1 m or x > 0.1 m), the beamforming levels are close. The mixing between surfacic and volumic focus points tends to favour the volumic ones. The GFs obtained for surfacic points have higher amplitudes resulting in lower beamforming power estimation. This amplification due to the reflection on the wing is also the reason why the maximum of beamforming using estimated GFs is not at the position of the source. However, after deconvolution, the estimated GFs enable to get the exact levels of a source located near the leading edge whereas the flow-corrected free-field ones give an over-estimation of its amplitude. Thus, acoustic imaging based on the estimated GFs enables the estimation of the amplitude and location of the acoustic source. The free-field GF with correction of the flow allows to obtain accurate results and also is a good approximation for doing imaging in this configuration.

Below the profile

We now focus on the plane 0.1 m below the wing which is discretized into 41×41 = 1681 focus points (blue area in figure 4.16). Beamforming results for a source positionned We also notice that the beamforming lobes downstream the position of the sources have higher amplitudes than those located upstream and (e)). This is because the trailing edge has a higher masking effect than the leading edge. In the expression of the conventional beamforming operator (1.18), the GF appears with a power of 2 in the numerator and a power of 4 in the denominator favoring masked regions.

This effect is even more important when the source is located in the visible region as shown in figure 4.22. For that case, secondary lobes have an amplitude greater than the level of the source and (e)). The correct level of 0 dB is found at the location of the source but the beamforming does not pass by a maximum at this position.

Discrepancies on the phases are not large enough to make the beamforming decrease because of the too important ratios between the GFs modulus. After deconvolution, the position and level of the source are found but artefact sources appear below the profile where beamforming levels are high (Figures 4.22 In every case, the total power of the source is 0 dB. In the figures 4.23-(b), the secondary lobes are still important but have been significantly decreased in figure 4.23-(c). Note that this solution is not always feasible as it depends on the characteristic sources to localize that can be narrow-band. Thus, it is better to avoid doing beamforming in volumic regions located not only in masked regions but also near flat reflecting surfaces.

The reflections on the surface of the wing do not have only negative effects here. Hence, the presence of the wing enables to increase the power of detection of the array in the Z-axis direction. If the focalization area is far enough from the surface to avoid this interference effect, the presence of reflecting surfaces can be an advantage. A possible application would be, in the framework of wind tunnel experiments, to intentionally place reflecting surfaces in order to improve the power of detection of a sensor array. However, those reflecting surfaces have to be positioned in a way that they do not perturbate the flow in the wind tunnel. For instance, it was numerically shown by [START_REF] Mimani | A simulation-based analysis of the effect of a reflecting surface on aeroacoustic time-reversal source characterization and comparison with beamforming[END_REF] that mounting a rigid surface on the floor of an anechoic wind tunnel can improve the capability of a line array of microphones to characterize experimental flow-induced noise sources.

The method also allowed to reduce the computational cost. If the GFs are required every 125 Hz, one period T = 8 ms of signal has to be recorded after the transient regime of t trans 3 ms. Using the monosource approach, one simulation of time t trans + T In this chapter, the proposed methodology was applied to two test cases representative for extern shapes for fuselage and wings: the diffraction by a rigid sphere and by a NACA0012 wing profile. For both cases, it enabled to get an accurate estimation of the acoustic GFs. When the number of focus points is greater than the number of microphones, we saw that the use of the reverse-flow reciprocity principle enables a decrease of the complexity of the estimation problem. The difficulty of the resolution of the resulting GF estimation problem is then driven by the number of microphone, usually much smaller than the number of focus points. Moreover, the reciprocity enabled to handle focus points located on the surface of diffraction objects because only pressure perturbations must be recorded at the focus points location. • A 9 x 9 m (29.5 x 29.5 ft) inlet featuring a dust filter, acoustic baffles, antiturbulence screens and a honeycomb. This part of the circuit aims to filter and straighten the flow in order to reproduce the flight conditions.

(a) (b) (c) ! "#$ %&'#(% !"# !$%"# !&# '%"# "# $%"# &# ($%"# ("# ('%"# (!&# (!$%"# (!"# (d) (e) (f) ! "#$ %&'#(% !"# !$%"# !&# '%"# "# $%"# &# ($%"# ("# ('%"# (!&# (!$%"# (!"#
• A contraction with an inlet diameter of 9 m (29.5 ft) and a nozzle exit diameter of 3 m or 2 m (9.8 or 6.6 ft). This part is in charge of the acceleration of the flow.

• An anechoic chamber (test section) in the shape of roughly a quarter of a sphere, with an internal radius of 9.6 m (31.5 ft). Walls are covered with absorbing foam which provide a very good anechoicity (minimised noise reflections), and low background noise in the 200 Hz to 80 kHz frequency bandwidth.

• A flow collector that aims to collect the flow after passing through the wind tunnel model. During experiments, one must ensure that the air jet coming from the inlet nozzle lands in the collector. Otherwise, it will lead to a consequent increase in noise generation and a flow recirculation.

• A diffuser.

• A fan silencer to avoid the fan noise from going upstream the flow and influence the measurement results.

• A centrifugal fan driven by a 7 MW asynchronous electric motor putting air into motion.

In the following, considered data were obtained in the framework of the TIMPAN cam- of the wing body is in the jet axis. The whole device is designed in order to minimize the noise generation by the additional structure. In the chosen configuration, the wing incidence of 11.5 • (13.5 • also available) is chosen. For this configuration, the tonal slat noise of the wing has more amplitude. Moreover, the lower jet deflection enabled a better collection of the flow on the collector resulting in a decrease of the background noise.

The other factor that influenced the choice for this configuration is the availability of data.

Acoustic instrumentation and measurement

For acoustic measurements, the wind tunnel is equipped with several microphones. One

!" #"

$"%"&"'("

)"%"*+"'(" ,"%"-"+./"(" ,"%"-"+.&"(" ,"%"+.&"(" ,"%"+./"(" 0"%"+./"(" 0"%"+.&"(" 0"%"-"+.&"(" 0"%"-"+./"("

(1'23453678" The first thing to do is to study the background noise and the one generated by mounting objects. These two types of noise must be negligible and if not, it is necessary to determine in which frequencies they are present. To ensure that, spectrum measurements with one microphone are performed in 3 configurations:

1. In the absence of the wing, only the mounting objects are presents.

2. In the presence of the wing, the slat is removed.

3. In the presence of the wing with the slat deployed.

The resulting spectra are shown in the figure 5.4. It appears that the noise measured by a single microphone is dominated by slat noise on most of the frequency range of interest except near 5.3 kHz. Around this frequency, spurious noise sources due to the support are very loud. Last but not least, it was ensured that the noise generated by the presence of the cross-shaped array is negligible.

Numerical simulation

CFD simulation

In the framework of the TIMPAN project, a 3D RANS computation of the flow around the wing and support was performed including the whole wind tunnel flow conditions The flow velocity magnitude in the XY-plane is represented in the figure 5.5. This plane is at position z = 0 m and contains the X-axis branch of the sensor array. At the exit of the nozzle, the flow is almost uniformly equal to 60 m/s. This velocity profile was verified by performing a comparison with experimental hot wire and Pitot measurements.

The diffusion of the velocity in the mixing layer during the propagation in the anechoic chamber is visible. An increase in the velocity in the extrados is visible. In this region, the velocity reaches 80 m/s. The flow recirculation in the slat cove and the thin shear layer is also visible. This shear layer is very important as it may causes instabilities during the CAA simulation step. As it will be seen in the next section, it is essential to filter the flow velocity field to use it as mean flow.

The flow velocity magnitude in the YZ-plane is represented in the figure 5.6. This plane is at position x = 1.932 m and contains the Z-axis branch of the sensor array The plane also contains the leading edge of the wing body and cuts the slat cove recirculation. It

shows that the flow velocity, circular at the output of the nozzle, has diffused. This validates that the whole support lies inside the velocity region. The low velocity behind the slat cove is due to the recirculation of the flow. To finish, this also shows that the sensor array is effectively outside the jet.

!"#$

!"#$% !"#% !"&'% !"&&% !"()% !"($% !"(% !"!'% !"!&% %&'()*+&,$-'+.$ /0.12$ 34).$

(a) CFL on the valid domain.

!"#$

!"#$% !"#% !"&'% !"&&% !"()% !"($% !"(% !"!'% !"!&% (b) CFL on the whole mesh.

CAA simulation

The CAA simulations were performed with ONERA's sAbrinA-v0 code on a body-fitted mesh. The useful part of the mesh corresponds to the region with -0.9 m ≤ x ≤ 1 m, -2 m ≤ y ≤ 0.2 m and -0.84 m ≤ z ≤ 0.84 m. The origin of the spatial system is at the middle leading edge of the wing body. In the region of uniform mesh, the spatial step is 4 mm. It is small enough to ensure more than 10 points per wavelength for the 1. The flow resulting from the CAA simulation is first interpolated on a coarse cartesian grid with a spatial step of 10 cm.

2. The flow on the coarse grid is then filtered by performing a convolution with 3×3 matrix with all the coefficients equal to 1/27. The filtering process is repeated another time in order to increase the filtering.

3. The flow on the coarse grid is then interpolated on the fine grid using a linear interpolation to obtain the flow of the figure 5.8-(b).

This filtering process has the advantage of being very easy to implement. If more filtering is needed, the step 2 can be repeated to increase the filtering. On the other hand, if less filtering is needed, the spatial step of the coarse grid in step 1 must be taken smaller.

However, the main drawback of this filtering technique is that it doesn't respect the condition of tangential flow on the surface. Considering a better filtering technique that respects boundary conditions may improve the quality of the CAA simulation. But, because of a lack of time and due to the important number of points in the mesh this method of filtering was privileged. The model of the flow obtained with this procedure is more complex than the cylindrical shear layer model flow shown in the figure 5.8-(c).

To be able to apply acoustic imaging technique, a focalization area must be chosen and discretized into focus points. The three regions of the focalization area and their discretization into N S = 12833 focus points is represented in the figure 5.9. The first region, colored in blue, is the intrados part of the airfoil and the slat (3150 focus points). The second region, colored in gray, extends the previous one by two plane surface upstream and downstream (respectively 3475 and 1650 focus points). A part of the support, representing the last region and colored in red, was also discretized (2×2279 focus points).

The procedure followed in this section is described in the figure 5.10. Acoustic GFs are estimated using a CAA simulation and are used to solve a problem of imaging with the

Acoustic imaging

Results

In this section, we present acoustic imaging results using CAA estimated Green's functions. The obtained acoustic maps will be compared with to the ones obtained with a cylindrical shear layer correction model. Hence, the theory on the effect of a thin

s i ⌅ s j (t) = 1 2 NF u=1 a (i,j) u cos ⇥ 2⇥f u t ⇤ (i) u + ⇤ (j) u ⇤ = 1 2 x (i,j) (t) (43)
i,n j,n 0 = (s i ⌅ s j ) (t n t n 0 )

⌅ i, j ⌅ J1, N S K n, n ⌅ J0, N V 1K c i,n = (s i ⌅ p M ) (t n ) ⌅ i ⌅ J1, N S K n ⌅ J0, N V 1K (44)
FFT {G M,i } (f ) =? (45) shear layer on the propagation of the sound has been widely studied in the literature. [START_REF] Amiet | Refraction of sound by a shear layer[END_REF] proposed a correction of the refraction effect of sound by shear layer. The method was validated afterwards by [START_REF] Bahr | Shear layer correction validation using a non-intrusive acoustic point source[END_REF]. To do so, they used a nonintrusive acoustic point source (pulsed laser system is used to generate a plasma point source) and showed that the method proposed by [START_REF] Amiet | Refraction of sound by a shear layer[END_REF] enables the recovering of the position of the source at several Mach numbers.

a u = C te (
In this section, the correction proposed by [START_REF] Elias | Correction des effets de réfraction dans la soufflerie sourde CEPRA19[END_REF] for the correction of the refraction effects in the CEPRA19 anechoic wind tunnel is used. The study of this correction procedure showed that it enables to recover precisely jet noise sources. In the following, beamforming is performed with the free-field Green's functions (1.14), the corrected Green's functions using the method proposed by [START_REF] Elias | Correction des effets de réfraction dans la soufflerie sourde CEPRA19[END_REF], and the CAA estimated Green's functions using the method presented in Chapter 2. A 3D view of the wind tunnel model, the microphone array, and the nozzle is given in the figure 5.12. In the following, the beamforming maps corresponds to the direction y < 0. The comparisons of the acoustic maps at 2 kHz and 3.8 kHz are presented in the figure 5.13. For these frequencies, the noise is due to the presence of the slat as shown in the figure 5.4. Hence, the measured power spectrum density is increased by more than 20 dB (at 2 kHz) and 10 dB (at 3.8 kHz) when the slat is present.

As expected, the use of the free-field Green's function is not adapted (cf. figure 5.13-(a) and (b)). The main lobe is larger at the frequency 2 kHz than at 3.8 kHz. This is in agreement with the theory that the main lobe varies as 1 /f bf with f bf is the considered frequency for beamforming. An offset of around d = 18 cm is present in the position of the detected acoustic sources. This is in agreement with the theory that gives d offset ≈ R nozzle M 1 + M 2 = 17.4 cm for the offset distance in the downstream direction where M = 0.176 is the Mach number. The corrected Green's functions using the cylindrical shear layer correction is shown in the figure 5.13-(c) and (d). For both frequencies, the generated noise is found to come from the slat cove.

The acoustic imaging results obtained using CAA estimated Green's functions are represented in the figures 5.13-(e) and (f). At 2 kHz, The noise source levels are slightly overestimated upstream. This is due to the mixing between focus points located on the surface of the wing and in the interior of the fluid. The overestimation corresponds to the main beamforming lobe and can also be seen at this frequency with corrected free-field Green's functions (figure 5.13-(c)). On the surface of the slat and the wing body, the size of the beamforming lobes and their amplitude differs completely. The lobes are more numerous but have a higher amplitude when compared to the free-field case. The width of beamforming lobes seems to depend on the curvature of the solid surface on which the focalization is done. Monopoles located on the surface see their resulting wavefront in the area of the microphone array impacted. However, we can say that globally, the estimated Green's functions are able to correct the effect of the flow.

Especially at 3.8 kHz, for this higher frequency, the beamforming lobes have smaller width and the source is found to be located in the slat.

In the following, the results based on the free-field Green's functions are not presented anymore. Acoustic maps obtained at 5.21 kHz, 5.39 kHz and 5.55 kHz are presented in the figure 5.14. In this figure, the results obtained by the corrected free-field Green's and (f)) are compared. As it was seen during the study of the noise spectra (Figure 5.4), at those frequencies, noise is present even when the wing and slat are not installed. Thus, the generated noise is likely to be coming from the support holding the wing. This is confirmed by the acoustic maps at 5.21 kHz and 5.39 kHz. At these frequencies, the main part of the noise is generated from the support. The slightly higher frequency 5.55 kHz corresponds to the transition from support noise to slat noise. Hence, for this frequency, the acoustic maps (Figure 5.14-(e) and (f)) show a mix from noise coming from the slat and the z < 0 part of the support.

At higher frequencies, the noise is mainly due to the slat noise. It is confirmed by the acoustic maps obtained at 6.07 kHz, 7.35 kHz and 8 kHz presented in the figure 5.15.

Discussion

The corrected Green's functions is well suited for the acoustic imaging configuration studied in this chapter. These Green's functions have been used for years in wind tunnel tests (refs) and are still in used today even on fan noise imaging [START_REF] Sijtsma | Using phased array beamforming to locate broadband noise sources inside a turbofan engine[END_REF].

For all the acoustic maps seen in the previous subsections, results obtained with CAA estimated Green's functions are consistent with those obtain using corrected Green's functions. However, results can still be improved.

On one hand, the maps obtained with the proposed method are noisy due to errors in the estimation of the Green's function. These errors may have their origins in several places.

More frequencies than N F = 601 must give more precise Green's functions and will result in less noisy acoustic maps. In the present case, acoustic paths from microphones to focus points are more numerous. The Green's functions are not as sparse as before. A longer time of simulation is necessary for the estimation of the Green's functions because more parameters are to be estimated.

On the second hand, several approximations made in the procedure of acoustic imaging. The geometry of the support was approximated by two parallelepipeds in order to facilitate both the mesh generation and the simulation process. Approximation of the geometry may have had an influence on the results. The approximation of the flow is also a source of error and must be as close as possible from the reality. However, it was necessary to filter the flow here in order to remove shear and boundary layers that create instabilities of the CAA simulation. The filtering is too blunt and the resulting flow doesn't respect the physics (mean flow penetrating the surface of the wing). A better filtering method, respecting boundary conditions is necessary to reduce the errors.

In this chapter, the proposed methodology was applied to the case of a high lift wing in the CEPRA19 anechoic wind tunnel. It constitutes the first case of validation of the method on industrial acoustic imaging configurations. The proposed method enabled to get an accurate localization of acoustic noise sources in the presence of a non potential mean flow (presence of a shear layer).

Conclusion

A methodology for the estimation of the Green's functions (GFs) based on numerical simulations for extern aeroacoustic imaging applications have been developed in this thesis. The method takes advantage of the sparsity of the GFs to minimize the computational cost and provide Green's functions that are sufficiently accurate to be used on realistic industrial configurations. It represents a continuation of Pene ( 2015)'s dissertation and various points were addressed.

The study of 3D cases has shown that the constraint on the l 2 -norm of the solution can lead to a deterioration of the solution. Two algorithms remain relevant for the resolution of the estimation problem. The lasso algorithm exhibits a behavior that is very stable and provides accurate GFs but suffers from a very slow convergence. On the other hand, the OMP algorithm provides GFs that are less accurate but its convergence is considerably faster and was privileged in this study. An improvement of the algorithms to avoid matrix inversion was also added for the acceleration of the algorithms. At each iteration, the inverse at the previous step is used to build the inverse of the augmented matrix. The addition of a cross-validation stopping criterion, more straightforward, also brought a consolidation of the method for avoiding overfitting.

A minimization algorithm for the generation of multisine signals with low pairwise crosscorrelation has been proposed. The algorithm allows the user to define the auto-power spectrum of the sources and iterates on the phases in order to minimize the common The proposed methodology was applied to test cases representative for extern shapes for fuselage and wings: the diffraction by a rigid sphere and by a NACA0012 wing profile.

Figure 1 :

 1 Figure 1: Measurements in the ONERA CEPRA19 anechoic wind tunnel.

Figure 2 :

 2 Figure 2: (left) An A340 flying over an array of microphones. (right) Beamforming maps.

  y∈V Θ(y, •) * g (y,0) (x, •) (t) d 3 y. (1.3) To obtain our model of the propagation of sound, we suppose that N S source points, located at positions (y 1 , • • • , y N S ), are excited as point monopole sources with signals (s 1 (t), • • • , s N S (t)). Only the pressure fluctuations at the microphone positions (x 1 , • • • , x N M ) are of interest. Then the forcing term reads as Θ(x, t) = N S i=1 δ (xy i ) s i (t), (1.4) and the pressure fluctuations at the location x m of microphone m is then p(x m , t) = N S i=1

  Figure 1.1: The propagation of sound is modeled by a LTI MISO system.

  φ(x) d 3 x.

Figure 1 . 2 :

 12 Figure 1.2: Acoustic imaging in the case of a source (a) in the Fraunhofer region or (b) in the Fresnel region of an array.

Figure 1 . 3 :

 13 Figure 1.3: (a) Angular beamforming with a rectangular window at 1, 2 and 4 kHz. (b) Angular beamforming at 1 kHz with a rectangular window, a Bartlett window and a Hanning window.

  ray pattern and consequently on the selectivity of the beamformer. The beamforming patterns at 1 kHz for a rectangular window, a Bartlett window and a Hanning window is presented in the figure 1.3-(b). The use of a window function to weight the receiver signals results enables a trade-off between the width of the main lobe and the height of the side lobes. A good overview of the techniques and methods for processing radar signals received by an array of antenna elements is the book of Haykin et al. (1993).

  b) CMM using free-field GFs.

Figure 1 . 4 :

 14 Figure 1.4: Focused beamforming in a homegenous infinite medium at rest. The 51 microphones are placed in the segment -0.5 m ≤ x m ≤ 0.5 m and y m = -0.5 m every 0.02 m. The source is placed at x s = (0, 0.5 m, 0) and its amplitude is 0 dB at the considered frequency f = 5 kHz. The GF of the medium is the free-field GF.

  d) CMM using analytical GFs.

Figure 1 . 5 :

 15 Figure 1.5: Focused beamforming in the presence of a shear layer. The 51 microphones are placed in the segment -0.5 m ≤ x m ≤ 0.5 m and y m = -0.5 m every 0.02 m.The source is placed at x s = (0, 0.5 m, 0) and its amplitude is 0 dB at the considered frequency f = 5 kHz. The medium is at rest in the bottom region (y<0) and the flow velocity is 100 m/s in the top region (y≥0). Acoustic maps obtained using CBF and CMM are shown using the free-field GF and a more precise GF that take into account the presence of the shear layer.

  d) CMM using analytical GFs.

Figure 1 . 6 :

 16 Figure1.6: Focused beamforming in the presence of a rigid sphere in a medium at rest. The 41 microphones are placed in the segment -0.4 m ≤ x m ≤ 0.4 m and y m = -0.4 m every 0.02 m. The source is placed at x s = (0.2 m, 0.3 m, 0) and its amplitude is 0 dB at the considered frequency f = 5 kHz. Acoustic maps obtained using CBF and CMM are shown using the free-field GF and the analytical GF that take into account the presence of the sphere.

  introduced a coordinate transformation by which the governing acoustic equations in a steady low-Mach number potential flow are transformed into an ordinary wave equation. The effect of a background flow field on sound propagation can be obtained by solving, independently, the potential (Laplace) equation for steady fluid flow, and the Helmholtz equation for acoustic-wave propagation in the transformed time variable, and then transforming the result into the physical time variable. This method was applied by Agarwal and Dowling (2007) in the study of the acoustic shielding by the silent aircraft airframe.

Figure 1 . 7 :

 17 Figure 1.7: Fan noise imaging on an installed turbofan. The propagation of acoustic waves is affected by the fan cowl and the flow making it difficult to localize acoustic sources.

  [START_REF] Candel | Résultats préliminaires sur la diffusion d'une onde acoustique par écoulement turbulent[END_REF] studied the scattering of spherical acoustic waves by a turbulent shear layer. For harmonic incident waves spectral analysis of the pressure field shows a fine structured tone broadening. Acoustic waves are strongly modulated in amplitude and phase during their migration in the turbulent mixing layer. Two lobes of scattered energy are found on both sides of the excited frequency. The study was led in the aim of characterizing the turbulent structures in the shear layer. The study of the space-time correlations and cross spectra between the direct signals and scattered fields allows an estimate of the convection velocity of the scattering vortices of the mixing layer. Koop et al. (2005) placed a test source within an open jet. The phase shift between the source and microphones outside the jet is measured. The dependency of the phase fluctuations on various parameters like frequency, flow velocity and microphone position is investigated. Kröber et al. (2010) addressed the problem of spatial coherence loss due to turbulent scattering and refraction. They proposed to use a calibration approach

  to compare acoustic results on a counter rotating open rotor CROR. Most dominant peaks in the spectrum of the rotor are well predicted by the method. The BEM thus doesn't meet the expectation to handle acoustic propagation through shear-layer encountered during open section wind tunnel tests.

Figure 1 . 8 :

 18 Figure 1.8: Procedure for Green's function estimation from a CAA simulation. During the simulation, focus points are excited using acoustic monopoles with excitation signals (s i ) i . During the simulation, pressure fluctuation signals (p caa m ) m at the microphones are recorded. The GFs are estimated by applying LTI system identification and adding a constraint on the solution. The estimated GFs can finally be used to perform acoustic beamforming.

  cover the frequency range of interest without numerical problems, excitations signals must be designed following several rules. In their survey of excitation signals for FFT based signal analyzers,[START_REF] Schoukens | Survey of excitation signals for FFT based signal analyzers[END_REF] explain that, to contribute constructively to the measurement, the energy must lie in the frequency range of interest. Moreover, the use of periodic signals measured over an integer number of periods enables to avoid leakage effects. If a non-integer number of periods are measured or if aperiodic signals are used, leakage effects cannot be avoided leading to a an increase of the needed measurement time to get a specified accuracy. In the following, multisine signals exciting this frequency range with harmonic components every ∆f are used:

  11) was proposed by Polifke et al. (2001) for N S = 2 non-periodic input signals for the reconstruction of acoustic transfer matrices in thermoacoustic systems. The methodology was also applied N S = 3 non-periodic signals by Föller et al. (2010) for the aeroacoustic characterization of T-junctions. It can be generalized for a larger number of periodic input signals. The Wiener filter Eq. (2the correlation matrix between input signals (s 1 , • • • , s N S ) and b is the crosscorrelation between input signals and pressure at the microphone: Γ in,jn

  .25) and is the most commonly used method of regularization of ill-posed problems. The hard constraint in (2.22) is often replaced by a soft constraint by using a Lagrange multiplier like in the Tikhonov (1943) regularization. The geometric representation of the ridge regression problem is shown in the figure 2.1-(a). The constraint region corresponds to

  non-convex. The geometric representation of the stepwise regression and l 0.8 -constrained problem are shown in figures 2.1-(b) and (c). For both of these cases, the solution can change completely with a small change on the covariates {u 1 , • • • , u N } and the observations v. The best solution estimate α can switch from one corner to the other.

  figure 2.3. In this test case, N S = 11 sources are located in the segment 0 ≤ x s ≤ 2 m at y=0 m every 0.2 m and z=0 m and the microphone is located at position (0, -2 m, 0).The OMP algorithm is used for the estimation of the GFs. Multisine signals (2.2) in the frequency range [4 kHz, 11 kHz] are used for the excitation of the sources. The procedure was applied for different number of frequencies N F = 141, 281 and 561 and the delimitation of the fitting (i) and overfitting (ii) phases relatively to the cross-validation stopping criterion has been highlighted. At each iteration of the algorithm, the GF estimates g m = g i,m i are updated in order to decrease the least squares criterion Q on the training set. Thus, the criterion Q decreases at every step making it difficult to differentiate fitting and overfitting phases. However, the criterion Q decreases, reaches minimum value, then increases when the algorithm starts overfitting. The best GF

Figure 2

 2 Figure 2.3: Cross-validation stopping criterion. The training (---) and validation (---) criteria are plotted as a function of the model complexity • 0 . Three numbers of frequencies N F =141, 281 and 561 are considered. The separation between the fitting and the overfitting phases has been highlighted.

Figure 2 . 4 :

 24 Figure 2.4: Time GFs obtained at (a) iteration 11 during the fitting phase, (b) at iteration 50 at the minimum of the validation criterion, and (c) at iteration 200 during the overfitting phase of the OMP algorithm. Vertical dashed lines represents the time of arrival of the first and the last echoes obtained analytically.

Figure 2

 2 Figure 2.5: Cross-validation stopping criterion. The curves were plotted for the OMP, the lasso, and the elastic net with λ 2 =0.001, 0.01, 0.1 and 1. The GFs between N S =11 focus points and a microphone are estimated. Multisine signals exciting N F =141 frequencies in the frequency range [4 kHz, 11 kHz] are used as excitation signals. Vertical dashed lines represent the separation between the fitting and the overfitting phases.

  very similar at iterations 200 and 300 (figure 2.6-(b) and (c)). The time domain GFs for the elastic net problem with λ 2 =1 are shown in the figure 2.7. Because of the constraint on the l 2 -norm, the algorithm gives solution that is less sparse. The amplitude of the coefficients is decreased because the energy is spread on more time components. This behaviour can be of interest when the Green's function is less sparse. In his PhD, Pene

Figure 2 . 6 :

 26 Figure 2.6: Time GFs obtained with the lasso algorithm at (a) iteration 11 during the fitting phase, (b) at the iteration 200, and (c) at iteration 300. Vertical dashed lines represents the time of arrival of the first and the last echoes obtained analytically.

Figure 2 . 7 :

 27 Figure 2.7: Time GFs obtained with the elastic-net algorithm at (a) iteration 11 during the fitting phase, (b) at iteration 189 at the minimum of the validation criterion, and (c) at iteration 300 during the overfitting phase of the elastic net algorithm. Vertical dashed lines represent the time of arrival of the first and the last echoes obtained analytically.

Figure 2 . 8 :

 28 Figure 2.8: Cross-validation stopping criterion. The curves were plotted for the OMP, the lasso, and the elastic net with λ 2 =0.001, 0.01, 0.1 and 1. The GFs between N S =51 focus points and a microphone are estimated. Multisine signals exciting N F =701 frequencies in the frequency range [4 kHz, 11 kHz] are used as excitation signals. Vertical dashed lines represents the separation between the fitting and the overfitting phases.

  et al. (1988) survey on excitation signals for FFT based signal analyzers. In the study of saturating systems, minimizing the crest-factor is of prior interest as it enables to use an input power as high as possible while keeping the device under test in its normal operation region. An algorithm for the CF minimization of multisine signals (3.1) was proposed by Guillaume et al. (1991). It is an iterative algorithm that gives a good estimate of the phases φ 1 , • • • , φ N F that minimize the Chebyshev norm:

  figure 3.1. Frequencies are chosen as f u = 1, 2, • • • , 40 Hz giving a period T = 1 s.Sources have a flat spectrum (a u,i = 1 ∀u, i). The algorithm is initialized with signals with random phases. The minimization is performed successively for p = 4, 8, 16, 32, 64 and 128. Only the phases of the multisine signal s 2 are optimized while the signal s 1 is kept unchanged. At the end of the process, the time aspect of signal s 2 is slightly modified and the crest-factor of the cross-correlation is decreased from 2.62 to 1.45. One can note that the energy (and thus the RMS value) of the cross-correlation given in (3.10) does not depend on the phases φ (i) u . Thus, the RMS value of the cross-correlation stays unchanged during the minimization process. When the number of frequencies is high, the signals vary very fast between ± RMS.

Figure 3

 3 Figure 3.2: Crest-factor evolution versus time at the application of the algorithm using a Levenberg-Marquardt algorithm (top) or gradient descent algorithm (bottom).Ten sources are considered with flat auto-power spectrum containing 100 frequencies f u = 1, • • • , 100 Hz. The algorithm is stopped when the crest-factor variation between two iterations becomes smaller than = 10 -4 two times in a row. These curves were computed on the same machine. Vertical dotted lines stand for 10 iterations of the algorithm and vertical dashed lines delimit the index p of the minimized norm l p .

Figure 3

 3 Figure 3.3: (Left) Crest-factors are plotted in function of the number of signals and the number of frequencies for random phases. Obtained values are between 3.06 and 6.25. (Right) Same figure for phases resulting of the application of the algorithm. Crestfactors have been decreased by more than a factor 2. Obtained values vary between 1.41 and 2.53.

  Fig. 3.3, the maximum of all the crest-factors obtained when doing the cross-correlation between two signals is plotted as a function of both the number of frequencies N F and the number of signals N S . Frequencies are chosen integers, that is to say f u = u ∀u ∈ 1, N F . The number of time steps is chosen as N T = 4N F , this is enough to satisfy the Shannon condition. The application of the algorithm enables to divide by a factor 2 the crest-factors in the considered range of values of N S .The obtained values are between 3.06 and 6.25 with random phases whereas they get between 1.41 and 2.53 after the optimization of phases. Moreover, for random phases, the maximum crest-factor is increasing with both N S and N F . After application of the minimization process, it appears that only the number of inputs N S really limits the minimization of the crestfactor.

Figure 3

 3 Figure3.4: (a) Time signal s 2 with random phases (---) and after (---) crosscorrelation minimization. Signal s 1 is unchanged so it is not represented here. (b) Crosscorrelation before (---) and after (---) application of the algorithm. The chosen region for low-correlation is represented in gray, here T M = 0.2 s. Cross-correlation crest-factor in the minimization zone is decreased from 2.62 to 0.12.

  .20) The minimization is done only on time steps in the interval [0, T M ] [T -T M , T [. Hence, it appears in (3.19) that minimization must be done for both positive and negative time delays. The part [T -T M , T [ comes from the fact that s

Figure 3 . 5 :

 35 Figure 3.5: Time aspects of a square (--), a sinus (--) and a triangle (--) signal. The three signals are T -periodic and have the same RMS value A.

  signal

  excitation signals. It appears in practice that, because of the fixed RMS value, the minimization of the Chebyshev norm also results in an augmentation of the l 1 -norm. To illustrate, the time aspects of three T -periodic signals are shown in the figure (3.5). The square signal s sqr , the sinus signal s sin and the triangle signal s tri have the same RMS value A and are expressed as:

4. 1

 1 Beamforming behind a rigid diffracting sphere 4.1.1 Experimental setup The data were obtained in the Aeroacoustic Windtunnel Braunschweig (AWB) anechoic wind tunnel (1.2 m × 0.8 m) of the German Aerospace Center (DLR). The experimental setup is presented in figure 4.1. The center of the 0.06 m radius sphere is taken as the

Figure 4

 4 Figure 4.1: (a): Diffraction by a rigid sphere in the AWB wind tunnel. (b): Schematic view.

Figure 4

 4 Figure 4.2: ONERA's SPARC acoustic source.

Figure 4 . 3 :

 43 Figure 4.3: Global mesh and CFL number.

Figure 4

 4 Figure 4.4: (a)-(b) components of the mean velocity field. (c) CFL number. Only the useful part of the mesh is shown meaning that the stretching zones are not plotted.

Figure 4 . 5 :

 45 Figure 4.5: Source and microphone signals.

Figure 4 . 6 :

 46 Figure 4.6: Spectra of one period of microphone signal (a) before and (b) after the convergence of the simulation.

Figure 4 . 7 :

 47 Figure 4.7: Diffraction of a point monopole source by a rigid sphere without flow (left) or in presence of a potential flow of 45 m/s (right). The source is represented with the violet diamond whereas microphone positions are represented with gray squares. The monopole source is excited using a multisine signal in the frequency range [4 kHz, 11 kHz] every ∆f=125 Hz (f exc 1 =4 kHz, f exc 2 =4.125 kHz, . . . , f exc 57 =11 kHz).

  figure 4.8-(a), (c) and (e) respectively for positions x m = -0.2, -0.07 and 0 m. The comparison is performed in terms of modulus and phase of the GFs in the frequency range [4 kHz, 11 kHz]. For the three microphone positions, a good agreement is obtained

Figure 4 . 8 :

 48 Figure 4.8: Comparison between GFs obtained experimentally (--), analytically (---) and using CAA simulation with one point monopole excited (+). Subfigures (a), (c) and (e) correspond respectively to microphone positions x m =-0.2, -0.07 and 0 m in absence of flow. Subfigures (b), (d) and (f) show the results obtained in presence of a mean flow of 45 m/s. For those last three subfigures, the GFs obtained by applying the proposed multisource strategy have also been represented.

  2 m]. These positions are highlighted in figure 4.10 where the pressure field resulting of the simultaneaous excitation of all the source points is presented. In order to obtain the training and validation sets S and S introduced in subsection 2.4, two simulations were performed. For the main simulation, corresponding to the training set S, exitation signals used are multisines as in Eq. (3.1) with N F = 1401 harmonics in the frequency range [f min , f max ] = [4 kHz, 11kHz] every ∆f = 5 Hz and with unitary auto-power spectrum:

Figure 4 .

 4 Figure 4.10: Presure field resulting by the simultaneous excitation of N S monopoles at the focus points (z s =+0.32 m). The perturbations of pressure are recorded at the microphones positions (z m =-0.3 m). The simulation is done in presence of a potential flow of 45 m/s. The sources are represented with the violet diamonds whereas microphone positions are represented with gray squares. The sources are excited using a multisine signals in the frequency range [4 kHz, 11 kHz] every ∆f = 5 Hz (f exc 1 =4 kHz, f exc 2 =4.005 kHz, . . . , f exc 1401 =11 kHz).

Figure

  Figure 4.11: Beamforming below a rigid sphere for three source positions (x s =-0.1, 0 and +0.1 m) and at three frequencies f=5, 7.5 and 10 kHz. The solid lines (--) corresponds to beamforming done with GFs estimated by the proposed method. Dashed (---) and dotted (• • • • • • ) lines correspond respectively to beamforming based on analytical GFs with and without taking into account the presence of the sphere. A correction was applied on analytical GFs in order to reproduce the effect of the 45 m/s potential flow.

Figure 4

 4 Figure 4.12: (a): Diffraction by a NACA0012 profile in the AWB wind tunnel. (b): Schematic view.

Figure 4 . 13 :

 413 Figure 4.13: Global mesh and CFL number.

  flow is first obtained by the use of the FreeFem++ software[START_REF] Hecht | New development in FreeFem++[END_REF]) because no analytical expression is availabe. The Laplace equation is solved on a large 2D rectangular domain in the X-Z directions meshed with triangles. The flow was then interpolated linearly on the CAA mesh and extended on the Y-direction using translational symmetry. The obtained potential flow is represented in the figure 4.14-(a) and (b). The mesh is cartesian outside the region with |x| ≤ 0.2 m and |z| ≤ 0.1 m containing the wing. A body-fitted mesh is realized in this region in order to apply hard wall BC on the surface of the wing (figure 4.14-(c)). The minimum mesh cell size is 1.4 mm and the time step is ∆t = 2.5 µs. The resulting maximum CFL number is 0.65. It is obtained near the surface at the maximum thickness of the wing and has to be smaller than 1 to ensure the stability.

Figure 4

 4 Figure 4.14: (a)-(b) components of the mean velocity field. (c) CFL number. Only the useful part of the mesh is shown meaning that the stretching zones are not plotted.

  are shown in figure 4.15. Three source and microphone positions are taken under consideration. The sources are placed in the visible upstream region (Figure 4.15-(b): x s = (-0.15 m, 0, -0.1 m) and x m = (-0.2 m, 0, +0.4 m)), exactly below the wing (Figure 4.15-(c): x s = (0, 0, -0.1 m) and x m = (-0.1 m, 0, +0.4 m)) and near the trailing edge of the wing (Figure 4.15-(d): x s = (+0.1 m, 0, -0.04 m) and x m = (0, 0, +0.4 m)).

  xs=(-0.15, 0, -0.1), xm=(-0.2, 0, +0.4), 55 m/s.

  xs=(0, 0, -0.1), xm=(-0.1, 0, +0.4), 55 m/s.

  xs=(+0.1, 0, -0.04), xm=(0, 0, +0.4), 55 m/s.

Figure 4 .

 4 Figure 4.15: Comparison between GFs obtained experimentally (--) and using CAA simulation with one point monopole excited (+) or by applying the proposed multisource strategy (•). Subfigures (b), (c) and (d) show the results obtained for three different positions in presence of a mean flow of 55 m/s.

  completely numerical test case because no experimental data are used. The setup is presented in figure 4.16-(a). The first surface, plotted in red, is in the plane z = 0 containing the wing and follows the extrados with |x| ≤ 0.2 m and |y| ≤ 0.2 m. The second surface, plotted in blue, is in the plane z = -0.1 m below the wing in the area |x| ≤ 0.2 m and |y| ≤ 0.2 m. The last surface, in green, is orthogonal to the wing profile at y = 0 with |x| ≤ 0.2 m and -0.2 m ≤ z ≤ 0.1 m. The three surfaces of focalization are discretized respectively into 1681, 1681 and 1291 points giving a total number of 4653 focus points. A cross-shaped sensor array consisting of 81 microphones is considered.

Figure 4 .

 4 Figure 4.16: (a) 3D representation of the three areas of interest and their discretization into 4653 focus points. The microphones of the cross-shaped microphone array used for acoustic imaging have also been plotted. (b) Pressure perturbation signal resulting from the excitation of 81 acoustic monopoles placed at the positions of the microphones of the array. The monopoles are excited using multisine signals in the frequency range [4 kHz, 11 kHz] every ∆f = 5 Hz (f exc 1 =4 kHz, f exc 2 =4.005 kHz, . . . , f exc 1401 =11 kHz). The mean flow of 55 m/s is reversed for the reciprocity principle to apply.

Figure 4 Figure 4

 44 Figure 4.17: Reverse-flow reciprocity principle. In the simulations, the monopole source is excited using a multisine signal in the frequency range [4 kHz, 11 kHz] every ∆f = 125 Hz (f exc 1 =4 kHz, f exc 2 =4.125 kHz, . . . , f exc 57 =11 kHz). (a) Direct approach. (b) Reverse approach. (c) Comparison of the GFs with experiments.

  .19. Acoustic GFs are estimated using a CAA simulation and are used to solve a problem of imaging with the classical beamforming method. The reverse-flow reciprocity principle is used in order to reduce the complexity of the GF estimation problem. During the simulation, microphones position are excited using acoustic monopoles with excitation signals (s m ) m and the pressure fluctuation signals (p caa i ) i at the focus points are recorded. Because no microphone array data is available, the aeroacoustic measurements are replaced by a CAA simulation to extract a synthetic cross-correlation matrix. The estimated GFs can finally be used to perform acoustic beamforming.

  Green's#func,on#es,ma,on#####(with#reciprocity)# resolu,on#of#the#GF#es,ma,on#problem######:# es,mated#Green's# func,ons# (p exp m ) m (b g i,m ) i,m Numerical#model# (with#reversed#flow)# # 4.#Acous,c#beamforming# acous,c#map#

Figure 4 . 19 :

 419 Figure 4.19: Procedure for acoustic beamforming based on CAA estimated GFs for the application on the NACA 2D wing. The reverse-flow reciprocity principle is used in order to reduce the complexity of the GF estimation problem. During the simulation, microphones position are excited using acoustic monopoles with excitation signals (s m ) m and the pressure fluctuation signals (p caai ) i at the focus points are recorded. The GFs are estimated by applying LTI system identification and adding a constraint on the l 0 -norm of the solution. Because no microphone array data is available, the aeroacoustic measurements are replaced by a CAA simulation to extract a synthetic cross-correlation matrix. The estimated GFs can finally be used to perform acoustic beamforming.

  simulation is shown in figure4.16-(b). In order to apply the reciprocity principle, the flow is reversed. Monopole sources are placed at the microphones positions and excited with multisinus (4.2) with N F = 1401 harmonics in the frequency range[4 kHz, 11 kHz] every 5 Hz. Pressure fluctuations are recorded at every focus point and an estimation problem of the form (4.4) is solved for each of them in order to extract the GFs. Three of the estimated GFs, corresponding to positions described in figure 4.15-(a), are plotted in figure 4.15-(b), (c) and (d). Again there is a good agreement with those obtained with the monosource strategy. The more defavorable case corresponds to figure 4.15-(b)when the source is centered below the wing. Hence, there is an important difference of the GFs amplitude between this focus point and the various microphones making problem (4.4) more difficult to solve.4.2.4.2 On the surface of the wingThe surface in the plane at z = 0 and following the extrados is discretized into 41 points in each direction giving a total of 1681 focus points (red area in figure4.16). The estimated GFs with the proposed multisource strategy are used to perform aeroacoustic imaging on the surface of a NACA0012 wing profile. Results of beamforming are presented in figure4.20 both in one and two dimensions. For the 1D-beamforming, only the streamwise brace of the sensor array is used in order to improve the results. Results after deconvolution have also been represented. The deconvolution was performed using the Least Angle Regression (LARS)[START_REF] Efron | Least angle regression[END_REF] model selection algorithm. The process is stopped when 99% of the energy of the acoustic map has been removed. It appears that the free-field GFs with flow is a good approximation. This is because a uniform flow is a good approximation in this case. In practice, the flow cannot always be considered uniform or potential and a mixing layer often appears between focus points and microphones. The proposed method however does not make the assumption of the flow being potential but only require a stationary flow.For a source near the leading edge of the wing, the use of the free-field GF leads to an error in the estimation of the position of the source. The source is estimated downstream its actual position due to the effect of the flow. The correction of the flow effect enables to localize the position of the acoustic source. It appears on the beamforming maps that the use of the free-field GFs with correction leads to a slight overestimation of the amplitude of the source. The amplification factor is due to the presence of a second acoustic ray diffracted by the leading-edge of the wing. The overestimation is more important at low frequencies for which the two rays reach the microphones with the same phase. On the surface of the wing (-0.1 m ≤ x ≤ 0.1 m), beamforming obtained with CAA-estimated GFs have levels 6 dB lower than those obtained with flow-corrected free-field GFs. The estimated GFs take into account the fact that, for focus points on the surface of the wing, two times the energy is scattered in the direction of the sensor array

Figure 4 .

 4 Figure 4.20: (a), (d): Acoustic beamforming results on the surface of the NACA0012 wing profile for a source positionned at (x s , y s , z s )=(-0.11 m, 0, 0) with unit power. The 2D acoustic maps are obtained using classical beamforming using CAA-estimated GFs. Positions of the source and microphones have also been represented. (b), (e): 1D beamforming results performed on the curve corresponding to y=0. Only the streamwise part of the array is used for 1D beamforming to improve the sensitivity. (c), (f): Results after deconvolution using the LARS algorithm. The considered frequency is 5 kHz ((a)-(c)) and 10 kHz ((d)-(f)).

Figure 4

 4 Figure 4.21: (a), (d): Acoustic beamforming results on a surface below the NACA0012 wing profile for a source positionned at (x s , y s , z s )=(0, 0, -0.1 m) with unit power. The 2D acoustic maps are obtained using classical beamforming using CAA-estimated GFs. Positions of the source and microphones have also been represented. (b), (e): 1D beamforming results performed on the curve corresponding to y=0. Only the streamwise part of the array is used for 1D beamforming to improve the sensitivity. (c), (f): Results after deconvolution using the LARS algorithm. The considered frequency is 5 kHz ((a)-(c)) and 10 kHz ((d)-(f)).

Figure 4 .

 4 Figure 4.22: (a), (d): Acoustic beamforming results below the NACA0012 wing profile for a source located at (x s , y s , z s )=(-0.15 m, 0, -0.1 m) with unit power. The 2D acoustic maps are obtained using classical beamforming using CAA-estimated GFs. Positions of the source and microphones have also been represented. (b), (e): 1D beamforming results performed on the curve corresponding to y=0. Only the streamwise part of the array is used for 1D beamforming to improve the sensitivity. (c), (f): Results after deconvolution using the LARS algorithm. The considered frequency is 5 kHz ((a)-(c)) and 10 kHz ((d)-(f)).

  -(c) and (f)). Free-field GFs with correction of the flow leads to an underestimation of the amplitude of the source and to the detection of a second source coming from the leading edge.4.2.4.4 Orthogonal to the wingThis time, we focus on the plane orthogonal to the wing which is discretized into 1291 focus points (green area in figure4.16). 2D beamforming results on the plane orthogonal to the wing are shown in figure 4.23. Only the streamwise part of the sensor array is used. With the flow-corrected free-field GF, an image source is detected due to the reflection of acoustic waves on the surface of the wing (Figure 4.23-(d)). Beamforming based on CAA-estimated GFs leads to very high amplitudes in regions where GFs have small amplitude. This happens in masked regions and near the surface of the wing (Figure 4.23-(a)). For a monopole placed near the surface, interferences between direct and image source lead to minimum and maximum of the modulus of the GF. However, the conventional beamforming operator (1.18) vary as the inverse of the square of the GF modulus. Performing beamforming near reflecting surfaces results in maximum and minimum of beamforming. The effect is more important near the surface as the direct and image sources have almost the same amplitudes leading to increasing interferences. This interference effect was studied by Sijtsma and Holthusen (2003) when studying the solution consisting in integrating image sources in the steering vectors in order to correct reflections during hard-wall closed-section wind tunnel experiments. It was also observed by Dougherty and Walker (2009) who performed broadband fan noise beamforming using steering vectors derived from annular hardwall duct modes. Broadband beamforming enables to decrease secondary lobes. When the frequency is changed, lobes do not appear at the same position whereas the position of the main lobe remains the same. Broadband beamforming was performed for a source emitting in third octave frequency range [4.47 kHz, 5.62 kHz] (Figures 4.23-(b) and (e)) and for a source emitting in the whole frequency range [4 kHz, 11 kHz] (Figures 4.23-(c) and (f)).

Figure 4 . 23 :

 423 Figure 4.23: Beamforming maps on the plane orthogonal to the NACA0012 wing profile for a source positionned at (x s , y s , z s )=(+0.1 m, 0, +0.1 m). Beamforming is done using estimated GFs ((a)-(c)) or free-field GFs with correction of the effect of the flow ((d)-(f)). Only the streamwise part of the sensor array is used to improve the sensitivity. Positions of the source and microphones have also been represented. The beamforming is done at the frequency 5 kHz ((a), (d)), in the octave band [4.47 kHz, 5.62 kHz] ((b), (e)) or in the whole frequency range [4 kHz, 11 kHz] ((c), (f)). In every case, the total power of the source is 0 dB.

Figure 5 . 1 :

 51 Figure 5.1: CEPRA19 open-jet anechoic wind tunnel.

  paign. It was interested in the experimental assessment of noise reduction concepts. A 2D high lift wing was placed at 2 m downstream the exit plane of the nozzle. Downstream the test measurement chamber the flow is run through a collector. The nozzle chosen in the experiments for the TIMPAN project is a 2 m diameter device allowing a maximum speed of 130 m/s.

Figure 5 . 2 :

 52 Figure 5.2: TIMPAN wind tunnel model.

  of the three main tasks of the EC TIMPAN project concerned High lift device activity. In particular, the both innovative concepts based on flow control technologies and mid-term noise reduction solutions as absorptive wing leading edge treatments were considered and a study of high-lift settings optimization through computational aero-acoustic methods was led. The microphones that are under interest in the present work are gathered in a crossshaped array of 41 microphone. The planar sensor array is used for acoustic source localization. It is placed in the plane y = -2 m. The geometry of the microphone array is shown in the figure 5.3. The cross-shaped antenna contains 41 half-inch condenser microphones numbered M1-M41. The positions of the various microphones are given in table 5.1. The central microphone faces the center of the leading edge of the wing body. During an acquisition session, 30 s of signal are recorded with all the microphones synchronized. They consist in 3.08 10 6 time samples collected at the sampling frequency F e = 102.4 kHz. Power and cross power spectral density estimates are then computed every 10 Hz in the frequency range [0 kHz, 20 kHz] using Welch (1967) overlapped segment averaging estimator furnished in the Matlab signal processing toolbox.

Figure 5 . 3 :

 53 Figure 5.3: Geometry of the microphone array.

Figure 5 . 4 :

 54 Figure 5.4: Pressure spectra in three configurations.

Figure 5 . 5 :

 55 Figure 5.5: Flow velocity magnitude in the XY plane.

Figure 5 . 6 :

 56 Figure 5.6: Flow velocity magnitude in the YZ plane.

Figure 5 . 7 :

 57 Figure 5.7: (a) CFL near the high lift wing profile. Only the useful part of the mesh is shown meaning that the stretching zones are not plotted. (b) Global mesh and CFL number.

Figure 5 . 8 :

 58 Figure 5.8: Surfaces iso-v based on the norm of the flow velocity are plotted at 30 m/s, 45 m/s and 60 m/s.

Figure 5 . 9 :

 59 Figure 5.9: Discretization of the focalization areas into 12833 focus points.

  Green's#func,on#es,ma,on#####(with#reciprocity)# resolu,on#of#the#GF#es,ma,on#problem######:# es,mated#Green's# func,ons# (p exp m ) m (b g i,m ) i,m

Figure 5 . 10 :

 510 Figure 5.10: Procedure for acoustic beamforming based on CAA estimated GFs for the slat noise imaging test case. The reverse-flow reciprocity principle is used in order to reduce the complexity of the GF estimation problem. During the simulation, microphones position are excited using acoustic monopoles with excitation signals (s m ) m and the pressure fluctuation signals (p caai ) i at the focus points are recorded. The GFs are estimated by applying LTI system identification and adding a constraint on the l 0 -norm of the solution. The estimated GFs can finally be used to perform acoustic beamforming.

Figure 5 .

 5 Figure 5.11: Pressure perturbation signal resulting from the excitation of 41 acoustic monopoles placed at the positions of the microphones of the sensor array. The monopoles are excited using multisine signals in the frequency range [2 kHz, 8 kHz] every ∆f = 10 Hz (f exc 1 =2 kHz, f exc 2 =2.01 kHz, . . . , f exc 601 =8 kHz). The mean flow is reversed for the reciprocity principle to apply.

Figure 5 .

 5 Figure 5.12: 3D view of the configuration.

Figure 5 . 13 :

 513 Figure 5.13: Beamforming maps at 2 kHz (left) and 3.8 kHz (right). The coherence with the monopole model operator (1.20) is plotted based on (a-b) the free-field GFs, (c-d) the flow-corrected free-field GFs, and (e-f) the CAA estimated GFs.

Figure 5 . 14 :

 514 Figure 5.14: Beamforming maps (a-b) at 5.21 kHz, (c-d) at 5.39 kHz and (e-f) at 5.55 kHz. The coherence with the monopole model operator (1.20) is plotted based on (left) the flow-corrected free-field GFs, and (right) the CAA estimated GFs.

Figure 5 . 15 :

 515 Figure 5.15: Beamforming maps (a-b) at 6.07 kHz, (c-d) at 7.35 kHz and (e-f) at 8 kHz. The coherence with the monopole model operator (1.20) is plotted based on (left) the flow-corrected free-field GFs, and (right) the CAA estimated GFs.

  l p -norm of the cross-correlation with increasing values of p. For low memory length systems, an adaptation of the algorithm to the generation of low-correlation zone (LCZ) signals has been presented. The results show that introducing an a priori information of the time interval on which cross-correlation has to be minimized enables to decrease it more effectively. When the reverberation time is known a priori, the condition number of the problem can be increased using LCZ signals because the cross-correlation can be decreased significantly. The benefits of using signals resulting from a minimization of the Chebyshev norm of the cross-correlation on the whole period require further investigations.

  46)(↵ g i,m ) i = argmin

					2
	(gi,m) i	p caa m	i	g i,m ⇥s i	2	s. t.	i	⇧g i,m ⇧ q ⇤ q	(47)
					2
	(↵ g i,m ) m = argmin (gi,m) m	p caa i	m	g i,m ⇥s m	2	s. t.	m	⇧g i,m ⇧ q ⇤ q	(48)
					2
	(↵ g i,m ) i = argmin (gi,m) i	p caa m	i	g i,m ⇥s i	2	s. t.	i	⇧g i,m ⇧ 0 ⇤ 0	(49)
					2
	(↵ g i,m ) m = argmin (gi,m) m	p caa i	m	g i,m ⇥s m	2	s. t.

m ⇧g i,m ⇧ 0 ⇤ 0 (50)

Acknowledgements xi

Introduction Acknowledgements

Chapter 5

Slat noise imaging using CAA-estimated Green's functions

In the previous chapter, we applied the Green's function estimation methodology on two test cases representative for extern shapes for fuselage and wings. By comparing the estimated Green's functions with analytical ones, we show that it enables to extract Green's function in the presence of both flow and diffracting surfaces. The method was also applied on synthesized microphone array data. Results showed that using high quality Green's functions can improve the quality of array data measurements.

In the present chapter, the method is applied in experimental microphone array data measurements. The details on the experimental setup are given in the section 5.1. Data were collected in the ONERA CEPRA19 anechoic wind tunnel [START_REF] Piccin | Cepra19: The ONERA large anechoic facility: A major tool for aeroacousitic measurements[END_REF]. The measurements were part of the European Commission (EC) Technologies to IMProve Airframe Noise (TIMPAN) project [START_REF] Manoha | Slat noise reduction using meshes: acoustic tests in CEPRA19 and numerical simulation of the 3D flow (EC project TIMPAN)[END_REF][START_REF] Perrin Decroux | TIMPAN -Technologies to IMProve Aiframe Noise[END_REF]. The numerical setup is given in section 5.2. The mean flow is obtained using CFD simulation solving the RANS equations and filtered afterwards (subsection 5.2.1).

It is then used during the CAA simulation in subsection 5.2.2. Again, the reciprocity will be used to decrease the required simulation time. In the section 5.3, beamforming maps using CAA estimated Green's functions are compared with ones obtained using classical shear layer correction method.

Experimental setup

Wind tunnel setup

The acoustic wind tunnel CEPRA19 [START_REF] Piccin | Cepra19: The ONERA large anechoic facility: A major tool for aeroacousitic measurements[END_REF] is operated by Onera in Saclay, France. It is an open circuit facility. The compressed air comes through a series of inlet devices:

Conclusion

For both cases, it enabled to get an accurate estimation of the acoustic GFs. In his perspectives, [START_REF] Pene | Antennerie numérique pour la caractérisation de sources aéroacoustiques en milieu complexe[END_REF] highlighted two conditions necessary for the application of the method in realistic industrial configurations: the method must be able to deal with a high number of focus (more than a thousand) and be able to deal with focus points located on walls. These issues were addressed by introducing the reverse-flow reciprocity: the acoustic response of a medium does not change when the source and microphone are swapped if the flow is inverted. As a result, the condition number of the problem is now driven by the number of microphones that rarely exceeds one hundred. The proposed methodology was then applied to the case of a high lift wing in the CEPRA19 anechoic wind tunnel. It is the first case of validation of the method on an industrial acoustic imaging configuration with around 13000 focus points. The proposed method enabled to get an accurate localization of acoustic noise sources in the presence of a non potential mean flow (presence of a shear layer).

For the geometry of the experimental envisaged in this work, the corrected free-field GFs is already a good approximation. The wind tunnel tests were not performed in the framework of this PhD and were designed in such a way that this condition is fulfilled. In order to validate the method, it must be applied on geometries of increasing complexity. This may require to conduct experiments specially dedicated to this purpose. The case of fan noise imaging still remains a challenge. For this case, the intrinsic nature of the acoustic GF is different and it may require to rethink the regularization strategy.

Appendix A

Validation of the free-field acoustic propagation

The 3D Free-field Green's functions in both time and frequency domains are:

where δ is the Dirac delta function and • the euclidean norm.

In Fig.

A.1, the pressure resulting from an harmonic source in a 3D free field is represented. Simulation was done using ONERA's sabrina v0 software. We obtain spherical wavefronts as the solution of this problem is the real part of the of the frequency Green's function at this frequency G (0,0) (x, f ). Appendix B

Acceleration of the cross-correlation minimization using FFT

This appendix presents the procedure to compute matrix

Computation of the matrix J T J

The element at line and column corresponding respectively to phases φ

v is:

Condition i = j means that only cross-correlation is minimized. Autocorrelation is not considered as it does not depend on the phases. The expression differs depending on if we are on a non-diagonal block k = l or on a diagonal one k = l.

excluding all the zero combinations. Using the notations:

3)

The expression becomes:

With the notations:

and denoting {z} the real part of a complex z, it comes:

where

u±v is the Discrete Fourier Transform (DFT) of the signal x (k,l) p-2 taken at frequency f u ± f v :

e -i(2π(fu±fv)tn) (B.8)

Thus the diagonal blocks of the matrix can be expressed in function of the non-diagonal ones:

Computation of the vector J T e

The element at line corresponding to phase φ (k)

u is:

Using notations (B.3) it comes:

where {z} stands for the imaginary part of a complex z. Again, if

the expression becomes:

where

is the DFT of the signal x (k,l) p-1 taken at frequency f u :

Hence, at each iteration of the algorithm, it is necessary to compute only one time the DFT X (k,l) u and X

(k,l) u of signals x (k,l) p-2 and x (k,l) p-1 respectively. The matrices can then be efficiently computed using equations (B.7), (B.11) and (B.15).

Appendix C

Diffraction of a monopole by a rigid sphere

In this appendix, we give the expression of the acoustic field diffracted by a monopole in the presence of rigid sphere. Suppose the acoustic source is at the position x i and the origin of the system of coordinates is the center of the sphere. Here, we are looking for the solution of the system of partial differential equations:

where G i,m (ω) = G i (x m , ω) is the frequency GF between the source at the position x i and the microphone at the position x m . The radius of the sphere is a and k 0 = ω/c 0 is the source wavenumber. With this convention, the free-field GF is expressed as:

However, taking into account the boundary condition in (C.1) modifies this solution.

The solution of (C.1) can be expressed as a sum of spherical harmonics:

where P l is the Legendre polynomial of degree l and θ i,m ∈ [0, π] is the angle between the acoustic source, the center of the sphere and the microphone. The amplitude G i,m (l, ω)

of the mode l is expressed in function of the spherical Bessel functions:

Appendix C Diffraction of a monopole by a rigid sphere where:

To determine the 4 constants in (C.5) we have to use continuity and boundary conditions.

To obtain only out-going waves in r m → +∞, A + i (k 0 , l) = 0. The expression of the last three constants is:

The result for a = 1, r i = 2 and k 0 = 4 where the sum in (C. ! Appendix D

Matrix inversion avoiding D.1 Active set increase

In this part, the case of an increase of the active set is considered. A new component

) is added to the active set A (q-1) = {(i 1 , k 1 ), . . . , (i N (q-1) , k N (q-1) )} giving the new active set A (q) .

The following relation can be used to compute the inverse of the new matrix using the previous one: q) invertible and Grammian Γ -1 A (q-1) exists and known

A (q-1)

is Grammian as: Γ A (q) (i,k),(j,l)

= (s i ⊗ s j ) (τ kτ l ) = ((s j ⊗ s i ) (τ lτ k )) * = Γ A (q) (j,l),(i,k) *

The vector v and the coefficient γ are defined as:      v = s i 1 ⊗ s i N (q) τ k 1τ k N (q) , . . . , s i N (q) -1 ⊗ s i N (q) τ k N (q) -1 τ k N (q) T γ = s i N (q) ⊗ s i N (q) (0)

Appendix D Matrix inversion avoiding

D.2 Active set decrease

On the other hand, if the active set is decreased, the following implication can be used in order to compute the inverse of the submatrix:

invertible and Grammian Γ

Often, the index to remove is not the last one added. In this case the formula is slightly different: