Y N Jaffré

T Aka-Ngnui

A Beroual

Non-Thermal Plasmas for NOx Treatment

Keywords: Corona Discharges, Energy, Non-Thermal Plasmas, Gas Treatments

This work is devoted to the determination of corona ignition threshold for non-thermal plasma generation and to the optimization of various kinds of plasma reactor geometries for exhaust gas treatment applications. The tested plasma reactor geometries were cylindrical. Some reactors have dielectric barrier made of glass or quartz in order to observe the discharges. First, the distributions of electric field and energy in reactors have been simulated. Then, current and voltage waveforms have been measured to check the discharge appearance. Geometrical singularities have been detected by both luminescence and current measurements to avoid unwanted current distortions. The experimental analyze shows the evolution of resistive and reactive components of current associated to corona discharges when increasing the voltage.

INTRODUCTION

Nitrogen oxides (NOx) are among the most important and critical pollutants that need to be reduced in exhaust gases. The present Selective Catalytic Reduction (SCR) processes can be improved for NOx reduction by a Non-Thermal Plasma (NTP) treatment. Plasmas are ionized gases, whose active species are generated by applying high electric fields in a plasma reactor. NTP alone does not reduce NO generated by a thermal engine to N2 and O2, but provides an oxidation of NO to NO 2 upstream of a traditional Selective Catalytic Reduction. The efficiency of this system depends on two main parameters: power supply and geometry. Plasmas have been studied for gas treatment since 1980. First applications were dedicated to power plants, for example in Riace, Italy, for ENEL's coal fired electrical generation plant (Civitano et al., 1986). Research for automotive transportation exhaust gas treatment started in the middle of the 1990's. Nowadays, the purpose is to assist the newest catalyst processes with non-thermal plasmas (NTP) to optimize the energetic consumptions and to improve the technology by mastering the discharge. This kind of treatment is of great importance since transportation exhaust gas regulations are becoming ever more stringent. European norm (standard specification) EURO 6 im-poses a reduction of 50% on automotive NOx emissions within 2015. This paper aims at the generation and optimization of plasmas for pollution control on automotive thermal engines. In the first part, we recall some fundamental features of discharge mechanisms. The second part lights on the influence of the electric field for NTP generation and link it to the electron energies. The third part shows results of external electric field simulations. And the last part presents respectively the experimental set-up and the measurements on NTP reactor.

DISCHARGE MECHANISMS

The current versus voltage characteristic is plotted on figure 1 for DC voltage supply with an interelectrode gap of 1 cm, plane to plane electrode geometry and standard air as dielectric. Three main regions can be observed and each region corresponds to a given mode. The first mode (dark zone) is characterized by the absence of luminous phenomena and by very weak currents. The applied voltage and thence the electric field induce elastic collisions and possible ionizations due to natural or artificial radiations. Since the discharge cannot be maintained by itself, it is called "non-autonomous". The second mode is characterized by a glow, a luminous activity defined on Figure 1 as the corona discharges. The latter starts when U exceeds U 0 ; thence the electric field brings sufficient energies to the electrons that ionize or dissociate molecules. Each electron newly generated contributes to electronic avalanches generation. Since the discharge propagates, the current increases, and the voltage falls down. Arc discharges involve high currents. The electrical energy consumption is the highest where arc discharges occur. The arc regime can be linked up with short circuits. Note that for the above conditions, the dielectric strength is of about 30kV/cm. Dark and Arc discharges have no interest for NTP generation. Dark discharges do not lead to oxidation or reduction processes. Arc discharges induce thermal mechanisms, high electrical energy consumptions and device destructions.

NTP PHENOMENA KNOWLEDGE

Our purpose is pollution control using corona discharges. These discharges are generated when the applied voltage (hence the electric field) exceeds a threshold value higher than the ionization and dissociation voltages of the considered gas. The electrons and ions The electrons move toward the anode and the ions (positive carriers) toward the cathode. By considering the NTP conditions, neutral molecules and ions are motionless compared to the electrons. The electron mobility can be approximated from the Chapman-Enskog theory of diffusion and simplified to:

µ = q mυ ( 1
)
where µ is the mobility, q the elementary charge (electron charge), m the mass of the electron and υ is the momentum transfer collision frequency. The mass of an electron is at least 2000 times less than any molecule. For an NTP, the external electric field mainly contributes to enhance electron energies. The electron drift velocities depend upon the magnitude, frequency and direction of the electric field:

ν = µ.E (2)
ν is the drift velocity and E is the external electric field.

The electrons are submitted to the Lorenz force which is:

F = q. E (3)
As concerns the current density j, it is given by the following relationship:

j = n e .q.µ.E (4)
where n e is the electron density. The kinetic energy of electrons is reliable to the enthalpy such as:

1 2 .m.ν = 3 2 .k.T e (5)
k is the Boltzmann constant and T e is the electron temperature expressed in eV (1eV=1.60.10 -19 J=11605K).

The statistical distribution of velocities of species is given by the Maxwell-Boltzmann statistical distribution equation. The electron energy distribution function can be defined as Maxwellian. Only a part of electrons can compete for dissociation and ionization processes. The equations above link the electron energy with the applied external electric field. By considering air at standard conditions of temperature and pressure, the physical mechanisms can be described by a series of complex equations (space charges, species interactions, transport, diffusion...). However, the external electric field remains a major parameter for NTP processes. Collisions occur according to parameters such as the energies, density of the species or collision cross sections. The mean free path of an electron is the drift distance between 2 electronic collisions:

λ = 1 n g .σ (6) λ = ν υ (7)
λ is the mean free path, n g is the molecular density of the surrounding gas and s is the cross section. For standard conditions in air, λ ≈1.10 -6 m; this can lead to either slowing electrons (elastic collisions) or generating electronic avalanche (inelastic collisions). The latter suggests that a satisfying electric field can produce adequate energetic electrons. More specifically, inelastic collisions lead to internal modifications of involved species. Some examples among numerous possible reactions are listed below:

e + O 2 → 2e + O + O + (8) e + O 2 → 2e + O + 2 (9) e + H 2 O → 2e + OH + H + (10) e + H 2 O → H -+ OH (11)
O -+ N O → N O 2 + e (12)
When considering a pure NO gas, the dissociation into N and O through inelastic collisions easily arises. The dissociation energy of NO is fixed at 6.5eV ([Pen93]); this value is obtained for a rough estimate electric field E=5.8MV.m -1 with λ=1.4.10 -6 m. The exhaust gases are mainly composed of N 2 and O 2 . Therefore, NTP does not generate a reduction but an oxidation. The oxidation takes place for energy levels above 5eV. The corona discharge ignition in the air can be calculated for cylindrical geometries from Peek formula

E c = 31.δ 1 + 0.308 √ δ.r (13)
with δ = 3, 92.p T g (14)
E c is the peak value of the corona field in (kV.cm -1 ), r is the radius of the conductor in (cm), p is the pressure in (cm.Hg) and Tg is the temperature in (Kelvin). The electric field distribution depends upon the applied voltage and the geometrical parameters. As the considered is coaxial, the electric field E can be calculated using Poisson's equation

∆.V = - q 0 (n i -n e -n n ) (15) E = -∇V (16) E(r, t) = U (t) r.ln Rout Rin ( 17
)
Where V is the voltage potential, n i and n n are ion and neutral densities, U is the differential applied potential, r is the distance parameter from internal electrode, R out is the external radius and R in is the internal radius. Both radicals O and OH are involved in oxidation processes of NO to NO 2 . These radicals are mainly produced by corona discharges in high electric field regions inducing energetic electrons according to the Maxwellian distribution of energy. But the oxidation itself takes place in low field regions during relaxing times of voltage applications or downstream of interelectrode gaps. In addition, some oxidations can appear upstream of interelectrode gaps resulting from discharge shockwaves.

CORONA DISCHARGE SIMULATIONS

Discharge on electrode geometries, magnitude and frequency of the voltage, electrical space charges, materials and gas parameters.

Geometry

In order to initiate the most efficient corona discharges for NTP generation, the geometry is highly divergent.

According to equation 16 it allows high electric field while voltage magnitudes remain below 25kV. The figure 3 shows a simulation of the internal radius influence on the electric field distribution for U=16kV and external radius of 15mm. Corona discharge ignition is calculated from equation 13. The curves show that a small radius enables to get a maximum value for E substantially higher than the corona ignition, but decreasing when increasing the gap from inner electrode. For a 15mm external radius, the internal radius is ranged between 0.2 and 4.2mm. While inner radius was set as variable on Figure 3, voltage U has been set as variable on Figure 3 where E was plotted for R in =0.35mm and R out =15mm. The corona ignition threshold is calculated for these geometric parameters: E corona =8.2kV/mm.

Power supply and arc limitation

The voltage has to be sustained below the disruptive voltage (U d ). However, the arc will not occur if the voltage application time over U d is too narrow. Streamers propagate with a non infinite velocity which means that the critical electric field can be reached by a pulse voltage during few nanoseconds but no arc. Care must be given on pulse repetition frequency, relaxing, rising and falling times. When relaxing time is too short and repetition frequency is too high, space charges can not be completely

Dielectric barrier

Another method to limit arc discharges consists of adding a dielectric barrier on one or both electrodes. Glass, quartz or any translucent materials are convenient for discharge observations as shown on figure 4. DC power supply cannot be used with Dielectric Barrier Discharges (DBD) because charges settle on barrier surfaces and cancel the electric field.

Materials and gas parameters

The NTP for NOx treatment has an oxidative leading function. Materials used for the electrodes have to remain unaffected by oxidation. Tungsten or stainless steel is preferred for building reactors. Exhaust gases can have rather high temperatures according to experimentations, with a maximum temperature of 550K. The insulating materials have to be matched to these temperature and corrosive requirements. The exhaust gas is simulated using an air compressor or generated using a 4strokes/1-cylinder fuel engine. Engine gas flow is ranged between 10l.min -1 and 150l.min -1 .

EXPERIMENTAL SETUP

The investigated system consists of a voltage source, a test reactor and an electrical measurement system. The voltage is measured thanks to a voltage divider, Tektronik 6015A (dividing ratio of 1000), with a bandwidth of 20MHz. The current is measured using 4 devices: 2 inductives, 1 resistive and 1 digital multimeter. The inductive devices are a transducer, Stangenes 0.5-0.1W, transforming current into voltage (ratio of 0.1V/A) with a bandwidth of 20MHz on one hand, a Rogowski coil associated to an integer (ratio of 0.1V/A) with a bandwidth range from 100kHz to 1GHz on the other hand. The resistive sensor is configured to obtain accurate measurement of the current on a 50Ω impedance. DMM device (Digital MultiMeter) is a NI PXI-4071 measuring current with ranges from 1pA up to 3A. This is mainly used to obtain the DC component and thus complete information on complex currents generated by discharges. In order to limit the current, a protective resistance is added to the circuit as shown on figure 5. Two types of voltage sources are used: (1) DC high voltage 60kV max; and

(2) AC high voltage 25kV (rms).

CORONA THRESHOLD MEASUREMENT

The measurements of high voltage currents are very sensitive. There are actual risks of arc discharges and thus the destruction of measurement devices. To avoid these shows almost a sin wave of current in phase with the voltage. Note that during experimentations we used a high speed CCD camera to observe discharges and eventual singularities (hot spot) that promote current peaks and arcs. The corona ignition threshold observation is in agreement with the previous calculation with an acceptable error of 5%. Luminance generated by the corona discharges appears for U=10.5kV (see figure 1). This would mean that singularities were avoided.

PROSPECTS AND CONCLUSIONS

The selective catalytic reduction (SCR) leads to NOx reduction. The products used for assisting the reactions are based on urea or hydrocarbons. The requirements to assist the reduction are the temperature and an oxidation stage converting a part of NO to NO 2 . The fastest reactions with NH 3 can take place for NO and NO 2 on stoichiometric ratios of 1:1. NH 3 is obtained by the urea hydrolysis. But, below 200 o C and during transient phases, the oxidative level doesn't work efficiently. The purpose is to create positive conditions for NOx reduction by SCR through the association with non-thermal plasmas whatever engine phases. This work contributes to discharge analysis for NTP generation. Corona discharges are well fitted for NTP reactors, but not for any voltage or geometry. In mastering electrical engineering, the combination of non-thermal plasmas and SCR has many additional advantages and promises to be a valuable economic way to treat pollutants in near future. Our researches lead to plasma technology experimentations with real exhaust gases using a thermal engine mounted on a test bench since we assumed that transient phases are not well transcripted by simulated gases.

Figure 1 :

 1 Figure 1: I-U characteristic for a DC discharge

Figure 2 :Figure 3 :

 23 Figure 2: E vs. inner radius for U=16kV and corona ignition value

Figure 4 :

 4 Figure 4: Corona Discharges with Glass (DBD)

Figure 5 :

 5 Figure 5: Equivalent Circuit of the Experimental Setup

Figure 6

 6

6 :

 6 Current and Voltage Characteristics of AC Corona Discharges for a Wire/Cylinder Geometry, No Dielectric Barrier. Umax ranges : (a) 9.5kV, (b) 11kV, (c) 14kV, (d) 16.5kV


N Bur
email: n.bur@estia.fr

P Joyot
email: p.joyot@estia.fr

On the Use of PGD for Optimal Control Applied to Automated Fibre Placement

Automated Fibre Placement (AFP) is an incipient manufacturing process for composite structures. Despite its conceptual simplicity it involves many complexities related to the necessity of melting the thermoplastic at the interface tape-substrate, ensuring the consolidation that needs the diffusion of molecules and control the residual stresses installation responsible of the residual deformations of the formed parts.

The optimisation of the process and the determination of the process window cannot be achieved in a traditional way since it requires a plethora of trials/errors or numerical simulations, because there are many parameters involved in the characterisation of the material and the process.

Using reduced order modelling such as the so called Proper Generalised Decomposition method, allows the construction of multi-parametric solution taking into account many parameters. This leads to virtual charts that can be explored on-line in real time in order to perform process optimisation or on-line simulation-based control. Thus, for a given set of parameters, determining the power leading to an optimal temperature becomes easy.

However, instead of controlling the power knowing the temperature field by particularizing an abacus, we propose here an approach based on optimal control: we solve by PGD a dual problem from heat equation and optimality criteria. To circumvent numerical issue due to ill-conditioned system, we propose an algorithm based on Uzawa's method. That way, we are able to solve the dual problem, setting the desired state as an extra-coordinate in the PGD framework. In a single computation, we get both the temperature field and the required heat flux to reach a parametric optimal temperature on a given zone.

INTRODUCTION

Automated Fibre Placement (AFP) is one the main technologies employed today to manufacture advanced composite laminates from unidirectional prepregs. This technique consists in laying and welding tapes of prepregs, building a laminate having the geometry and the desired mechanical characteristics.

This process has been widely studied [1,2,3,4] with numerical models becoming more accurate, robust and complex. New techniques of computation [5,6], based on variables' separation, have enabled models' enrichment by adding parameters as extra-coordinates [7]. Thus, the Proper Generalised Decomposition (PGD) method leads to multi-parametric virtual charts providing a whole set of solutions for each combination of the considered parameters [8,9,10,11,12].

Then, the computational vademecum can be exploited on-line for process control or process optimisation purposes. Indeed, within the AFP we want to efficiently control the heating power: tapes have to be heated enough to ensure the melting of the matrix coating the fibres and the cohesion with the previously laid tapes, while not exceeding a threshold from which material burns.

Therefore, we took advantage of the PGD to build off-line virtual charts in order to determine the best power associated to a draping velocity profile [13].

However, in those simulations, solutions were computed from equations provided by underlying physics of the studied phenomena, the optimisation being carried out as post-process.

We propose here to compute directly the solution of an optimisation problem in order to get the separated representations of both the field and the control to obtain it. That is to say the optimisation is made directly off-line, reducing the cost of the post-process and improving the real-time control of the AFP.

Within the next section, we present the equations governing the phenomenon under consideration. In order to have a reference solution the system is solved by standard finite element method (FEM). Thereafter, section b) focuses on the writing and solving of the optimal system. We improve the obtained results by applying Uzawa's method in section b). Lastly section b) addresses few conclusions and perspectives.

PROCESS MODELLING

AFP process can be modelled with an heat equation associated with the next boundary conditions. Heat source applies on boundary Γ P (see Figure 1); considering a wide enough domain Ω, we set a Dirichlet's condition taking into account a fixed temperature on Γ R ; then, for the sake of simplicity, an homogeneous Neumann's condition is applied on others boundaries.

Γ B Γ L Γ T Γ P Φ Ω FIGURE 1. Domain of study.
This leads to the following advection-diffusion equation 1

                 -div (K∇u) + ρC p V∇u = 0 in Ω; u = 0 on Γ R ; -K∂ n u = 0 on Γ B ∪ Γ L ∪ Γ T ; -K∂ n u = -Φ on Γ P ; (1)
with

K = k 0 0 k ⊥ et V = v 0 T .
This system can obviously be solved in a classical way with standard FEM. However, since we want a reference solution to latter compare with results from other methods, we have to implement comparable systems. Consequently, in order to compare solution in separated representation, we take advantage of the tensor product method introduced by R. E. Lynche, J. R. Rice and D. H. Thomas [14] [15]. The writing of discrete form of Equation ( 1) separating the two space coordinates and using tensor forms of operators allows the handling of shape functions in the separated representation, the system remaining solved with a global, non separated FEM.

We can thus obtain the temperature field from the input Φ modelling the laser heat flux. For some couples power/speed we compute the corresponding reference solutions.

From these results we extract the temperature profiles on boundary Γ P where the source is applied. Indeed, we want to control the temperature on this same boundary during the process, to reach an optimum, providing the ideal heat flux. We compute also P S ray , the power seen by Γ P , as the integrand of the flux on this boundary multiplied by the width of the shone surface S ray .

We gather in Table 1 some key values for these four reference solutions, for the purpose of latter comparison. Due to velocity values, we increase the number of nodes to contain the Péclet number, avoiding the requirement of a stabilisation technique

v = 0.001 m • s -1 v = 0.01 m • s -1 v = 0.1 m • s -1 v = 1 m • s -1 Pw = 600 W Pw = 1897 W Pw = 6000 W Pw =

SETTING UP THE OPTIMAL SYSTEM

As announced in the Introduction, The PGD method provides multi-parametric virtual charts that can be used to control the AFP process [13].

Another way to go on is to use the optimal control theory. Within the AFP process, a heat flux provided by a laser melt the thermoplastic. The difficulty consists in determining the best power of the laser to reach an optimal temperature to melt the thermoplastic enough but not too much.

Thus we consider the following cost-function to be minimised, since the flux is applied only on Γ P , part of the boundary

J(u, Φ) = 1 2 Γ P (u -u d ) 2 + αΦ 2 dγ, (2)
subject to the state advection-diffusion equation ( 1), with α cost parameter of the command. That way Φ is used as control and we want to reach u d on boundary Γ P . The domain of study is depicted on Figure 1.

The corresponding Lagrangian writes

L (u, Φ, p) = J(u, Φ) + Ω p -div (K∇u) + ρC p V∇u dx. ( 3
)
To find a stationary point of the Lagrangian we set

           ∇ u L (u, Φ, p) = 0; ∇ Φ L (u, Φ, p) = 0; ∇ p L (u, Φ, p) = 0. ( 4
)
Expanding these equations leads to the non-linear optimality system (see [START_REF] Lions | Optimal Control of Systems Governed by Partial Differential Equations[END_REF]), whom weak form writes, with test functions u et p ,

       Ω K∇u • ∇u + Ω ρC p V∇uu -Γ P 1 α pu dγ = 0 ∀u ; Ω K∇p • ∇p + Ω ρC p Vp∇p + Γ P up dγ = Γ P u d p dγ ∀p . (5)
In Equation ( 5), fields u and p are coupled. The computation of such a problem can be achieved using mixed formulation, within a standard FEM as well as in the PGD framework as described thereafter.

The discrete form of the state variable u (x, z) and the adjoint parameter p (x, z) are expressed in tensor product form

U = ∞ i U i x ⊗ U i z and P = ∞ i P i x ⊗ P i z . (6)
The vector Ψ which brings together nodal values of u and p takes the form

Ψ = U P = ∞ i U i x ⊗ U i z P i x ⊗ P i z . (7)
The discretised weak form of Equation ( 5) is written

Ψ T AΨ = Ψ T B were A = 8 i=1 A j
x ⊗ A j z and B = B x ⊗ B z will be expressed in a tensor form, and with the test function defined by

Ψ = U x ⊗ U z + U x ⊗ U z P x ⊗ P z + P x ⊗ P z . ( 8
)
This discretised problem can then be solved with the PGD method, for different values of the velocity, the desired state u d coming from the corresponding reference solution.

Since the goal of optimal control is to minimise the distance between the unknown u and the desired state u d on the boundary Γ P , we retrieve the wanted temperature profile on this boundary. However, to reach it, the flux, computed by PGD, to be applied take substantially away compared to the FEM input, as speed increases (see Figure 2).

v = 0.001 ms -1 v = 0.1 ms -1 v = 0.01 ms -1 v = 1 ms -1

UN-MIXING THE OPTIMALITY SYSTEM

Applying the PGD framework on the optimality system (5) can produce, in our case, weird results. Instead of searching a workaround within numerical stabilisation techniques, we propose to solve this system without mixed elements. For this purpose, we use Uzawa technique [17] [18].

Given an initial guess p (0) for p, Uzawa's method consists in our case of the following coupled iteration:

       A u u (k+1) = A pu p (k) p (k+1) = p (k) + ω A pu u (k+1) -u d + A p p (k) (9)
where ω > 0 is a relaxation parameter. Tensors are defined by 2 summarises the previous simulations, collecting some key-values. Thus, solving the optimality system with the PGD algorithm, but without mixed elements using Uzawa method leads to consistent results, since both temperature and heat flux remain similar to those expected (using standard FEM).

A u = k K u x ⊗ M u z + M u x ⊗ k ⊥ K u z + ρC p vH u x ⊗ M u z A pu = - 1 α M pu x ⊗ M pu

FIGURE 2 .

 2 FIGURE 2. Flux on Γ P for different speeds.

z

  A p = k K p x ⊗ M p z + M p x ⊗ k ⊥ K p z -ρC p vH p x ⊗ M p z A up = M up x ⊗ M up zThis coupled iteration is computed within a fixed-point loop, each field u and p being solved by PGD. As previously, Figure3shows the computed flux to reach the desired temperature u d on the boundary Γ P .v = 0.001 ms -1 v = 0.1 ms -1 v = 0.01 ms -1 v = 1 ms -1

FIGURE 3 .

 3 FIGURE 3. Flux on Γ P with Uzawa's method Table2summarises the previous simulations, collecting some key-values. Thus, solving the optimality system with the PGD algorithm, but without mixed elements using Uzawa method leads to consistent results, since both temperature and heat flux remain similar to those expected (using standard FEM).

TABLE 1 .

 1 Key values for reference solutions

x z l 1 l 2 l 3 h Γ R

CONCLUSION

In order to control the AFP process, we proposed a method based on the PGD technique allowing the inclusion of parameters as extra-coordinates. Due to numerical instabilities, we transformed a mixed formulation to a coupled problem, improving significantly the results.

Next step consists in considering the desired state u d as a new coordinate. Thus we will be able to build directly the process control as a multi-parametric virtual chart.


S Maugeais
email: sylvain.maugeais@univ-lemans.fr

J Gilbert

Brass player's mask parameters obtained by inverse method

An optimization method is proposed to find mask parameters of a brass player coming from a one degree of freedom lip model, with only constant mouth pressure and periodic mouthpiece pressure as input data, and a cost function relying on the waveform and the frequency of the signal. It delivers a set of parameters called C -admissible, which is a subset of all mask parameters that allow the inverse problem to be well defined up to an acceptable precision. Values for the mask parameters are found that give a good aproximation of real signals, with an error on the playing frequency of less than 5 cents for some notes. The evolution of the mask parameters is assessed during recordings with real musicians playing bend notes and their effects on the playing frequency are compared to the theoretical change on a model.

Introduction

Models of brass instruments have been used for a long time to understand the physics (cf. [START_REF] Campbell | The Science of Brass Instruments[END_REF]) and synthesize their sound (cf. [START_REF] Harrison-Harsley | Physical Modelling of Brass Instruments using Finite-Difference Time-Domain Methods[END_REF]). However, one remaining difficulty is the calibration of these models as many of the constants appearing in the mathematical equations are difficult to measure, in particular when they involve human body parts such as the lips of a brass player. A good calibration of these parameters is useful not only for sound synthesis, where a realistic sound is the main goal, but also for theoretical purposes, as the equations governing models of brass instruments are nonlinear, and are therefore very sensitive to changes in the constants of the model.

In the literature, many models exist for the lips, with different degrees of complexity, either as a one degree of freedom oscillator (which is the most common), or a two degrees of freedom oscillator taking into account different polarities (cf. [START_REF] Boutin | Trombone lip mechanics with inertive and compliant loads ("lipping up and down")[END_REF]), and models trying to come closer to the geometry of the opening section of vibrating lips (cf [START_REF] Campbell | The Science of Brass Instruments[END_REF], section 5.1.2). However, the more complex the model, the more constants there are that need to be calibrated, and the higher the uncertainty. The question of calibration of the parameters is not restricted to lip models of brass players, and a related problem is that of the embouchure of reed instruments which may seem easier as measures can be undertaken directly on the reed. The study of the reed parameters is the source of a large literature (see for example [START_REF] Avanzini | Modelling the Mechanical Response of the Reed-Mouthpiece-Lip System of a Clarinet. Part I. A One-Dimensional Distributed Model[END_REF][START_REF] Van Walstijn | Modelling the Mechanical Response of the Reed-Mouthpiece-Lip System of a Clarinet. Part II: A Lumped Model Approximation[END_REF][START_REF] Muñoz Arancón | Estimation of saxophone reed parameters during playing[END_REF][START_REF] Chatziioannou | Estimation of Clarinet Reed Parameters by Inverse Modelling[END_REF][START_REF] Helie | Inversion of a physical model of a trumpet[END_REF][START_REF] Chatziioannou | Investigating Clarinet Articulation Using a Physical Model and an Artificial Blowing Machine[END_REF][START_REF] Smyth | Toward an estimation of the clarinet reed pulse from instrument performance[END_REF]) which can be used as a source for methods dedicated to the lips.

Concerning brass instruments, the body of literature is more reduced, although many articles deserve to be cited and will serve as reference for the present work (see [START_REF] Velut | How Well Can Linear Stability Analysis Predict the Behaviour of an Outward-Striking Valve Brass Instrument Model?[END_REF] for a table summarizing known values). Most notably, [START_REF] Elliott | Regeneration in brass wind instruments[END_REF] who was one of the first to give a complete set of parameters, [START_REF] Brigitte D'andréa | Asymptotic State Observers for a Simplified Brass Instrument Model[END_REF] who built an asymptotic state observer, [START_REF] Vergez | The BRASS Project, from Physical Models to Virtual Musical Instruments: Playability Issues[END_REF] who used simulated annealing with a cost function depending on playing frequency only, [START_REF] Boutin | Trombone lip mechanics with inertive and compliant loads ("lipping up and down")[END_REF] with data coming from high speed camera and [START_REF] Fréour | Parameter of a physical model of brass instruments by constrained continuation[END_REF] using bifurcation diagrams. The present article aims to provide a new method to identify embouchure parameters that can be used with very little apparatus on actual musicians, and that can follow their evolution while playing. It introduces a new cost function which is a combination of [START_REF] Vergez | The BRASS Project, from Physical Models to Virtual Musical Instruments: Playability Issues[END_REF] (for the frequency part) and [START_REF] Chatziioannou | Investigating Clarinet Articulation Using a Physical Model and an Artificial Blowing Machine[END_REF] (without the displacement), together with some penalization (see section 2.4). An optimization algorithm is used to minimize this cost function on recordings with actual musicians and the results are discussed (cf. section 3).

Method

Model

To reduce the amount of parameters that must be calibrated, the model chosen in the present article is the simplest one described by (2.1), with only three equations (cf. [START_REF] Mattéoli | Minimal blowing pressure allowing periodic oscillations in a model of bass brass instruments[END_REF]), that proved to replicate many of the properties of brass instruments (see for example [START_REF] Mattéoli | Diversity of ghost notes in tubas, euphoniums and saxhorns[END_REF]). It relates the mouth pressure p m to the mouthpiece pressure p and the opening h through a spring-mass-dashpot equation describing the lips, a valve effect computing the flow u through the lip from the difference of pressure, and the expression of the input impedance:

                 ḧ + ω ℓ Q ℓ ḣ + ω 2 ℓ (h -H) = pm-p µ ṗn = Z c C n u + s n p n p = 2ℜ N n=1 p n u = wh + sgn(p m -p) 2|pm-p| ρ (2.1)
where h + = max(h, 0), w is the width of the lips, and Z c the characteristic impedance.

Here, the impedance is decomposed using modal analysis as a sum of simple fractions

Z(ω) = Z c n C n jω -s n + C * n jω -s * n (2.2)
and p n being complex valued. The variables p n form a decomposition of the mouthpiece pressure using the special form of the impedance. They are not really modal coordinates, as they are not naturally orthogonal for some scalar product, but give a convenient way to solve the problem. The input impedance is measured through an impedance bridge and is therefore known. It is decomposed into formula (2.2) using the Rational Fraction Polynomial method (see [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF], section 4.4.3), giving the values for s n and C n which are therefore fixed for a given instrument. The value of C n and s n computed for the Bb trombone in first position and used in the measures (basse trombone Courtois and mouthpiece Holton) are given in table 5. The mask parameters are the remaining constants appearing in equations (2.1), and are the lip angular resonance frequency ω ℓ = 2πF ℓ , the quality factor Q ℓ , its surface density µ and its opening at rest H. The set of equations 2.1 can be rewritten as a real valued ordinary differential equation (cf. [START_REF] Mattéoli | Minimal blowing pressure allowing periodic oscillations in a model of bass brass instruments[END_REF]) ẏ = f (y) of order 1 in dimension 2N + 2 with a state variable

y = [h, ḣ, ℜ(p 1 ), Im (p 1 ) . . . , ℜ(p N ), Im (p N )]. (2.3)
Such a formulation is very convenient both for time numerical simulation, and for the study of bifurcation diagrams, as computed by the software auto-07p [START_REF] Doedel | Auto-07p, continuation and bifurcation software for ordinary differential equations[END_REF] (see section 2.9). The time simulations performed in this project are based on the Runge-Kutta 4 algorithm and are coded in C language as a large number of them are performed. The sampling rate is set to 44100Hz as it gives a sufficient precision for the results (comparable to measured data), ensures the stability and convergence of the numerical scheme, and is sufficiently fast (0.18s for a one second signal). The time simulations have to be performed on a time range long enough so that stationary regime is attained. In practice, this means signals of up to 4 seconds have to be simulated as transient regime can be quite long for some mask parameters.

For a set of mask parameters

M = (ω ℓ , Q ℓ , H, µ), (2.4)
we write p M for the mouthpiece pressure obtained by solving the system (2.1) with zero initial value

y(0) = [0, • • • , 0] (cf. equation 2.
3). All the other parameters of the model, including p m , are fixed and constant during one optimization.

Experimental protocol

The goal of the project is to get access to as many mask parameters as possible with as little apparatus as necessary so as not to hinder the musician's playing. We therefore focused on only two piezoresistive pressure sensors (Endevco 8507C-5): one in the mouth of the musician (or artificial mouth), and one in the mouthpiece. Simultaneous recordings give access at each time step to the quantities p m and p. Both signals are sampled at 44100Hz. As only one period of p is needed by the algorithm below, the signal can be broken into pieces to see the evolution of parameters (see section 2.4). Once a period is chosen, p m is averaged on the same time period to have a constant value to feed the numerical simulation.

Three sets of measures were performed with two experienced amateur trombone players, one of them being recorded twice, labeled A1, A2 and B in the rest of this article. For each session, the musician was asked to play 6 notes on a Bb bass trombone in first position (Bb2, F3, Bb3, D4, F4, Bb4), together with a bend on F4 (first down, then up), and a crescendo on F4.

Extraction of a signal's period

For a given sampled signal p, either simulated (as in equation (2.1)) or measured (see section 2.2), the determination of a period is critical for the method and the first step to perform it is the identification of the periodic regime and its frequency. This is done using a python implementation [START_REF] Guyot | Fast python implementation of the yin algorithm[END_REF] of the Yin algorithm (see [START_REF] De | YIN, a fundamental frequency estimator for speech and music[END_REF]). The Yin algorithm produces an estimator called the harmonic rate, which is a real number between 0 and 1, that gives a quantification of how periodic the signal is: A harmonic rate very small (ideally 0) meaning that the signal is close to being periodic. In practice, we consider the signal to be periodic when the harmonic rate is lower than 10 -3 , and extract the part of the periodic regime with the smallest harmonic rate.

Yin also gives an estimation of the instantaneous frequency F at each point. From this, it is already possible to extract a waveform p † of duration exactly one period in the periodic regime. However, as this waveform has to be compared to another one coming from a reference signal, a phase condition has to be fixed. We therefore demand that all the waveforms • begin by crossing 0,

• in an increasing way.

(2.5) This is always possible in practice as p has a mean value equal to 0. The normalization is achieved by considering the waveform p † and shifting it to the left until the first point satisfies the phase condition, giving rise to a new waveform p which is used as a reference.

It should be noted that on a general signal, the phase condition may not be sufficient to uniquely determine p. However, for the signals obtained either numerically or experimentally, this condition proved to be sufficient.

For a set of mask parameters M, we also write F M and p M for the frequency and normalized waveform of the signal p M obtained in section 2.1.

Definition of the cost function

The goal of the cost function is to try to compare two periodic signals p ref and p, from which frequencies F ref , F and waveforms p ref , p are extracted. The signal p ref can be either a recorded signal, or a simulated signal obtained from known mask parameters (for test purposes), and is the reference against which the model outputs are compared. Although in theory it should be enough to compare p ref and p, it puts too much emphasis on the waveform itself, and too little on the frequencies. As both timbre and intonation are important for the applications, it is necessary to add an extra weight to the difference in frequencies. The preliminary cost function C • is therefore

C • (p ref , p) = 1 ∥ p ref ∥ 2 2 min 1 F ref , 1 F 0 ( p ref (t) -p(t)) 2 dt + α F (1200 log 2 (F ref /F )) 2 (2.6) with ∥ p ref ∥ 2 2 = 1/F ref 0 p ref (t) 2 dt.
The first term of the sum is the square of the relative RMS difference, and the second one the square of the relative frequency difference.

The choice of the constant α F guides the optimization procedure either toward a better approximation of the waveform (small α F ) or toward a better approximation of the frequency (big α F ).

In our case, the choice of α F = 0.02, obtained by trial and error, leads to good results during optimization for trombone sounds, in that intonation (errors around 10 cents) and waveforms (errors around 30%) are respected. In particular it means that the difference in cents between two signals is

≤ C • (p ref , p)/α F .

Penalization

As is already well known (cf. [START_REF] Helie | Inversion of a physical model of a trumpet[END_REF]) the inversion problem is not well defined and it is actually easy to find multiple sets of mask parameters which give signals with very similar waveforms (see table [START_REF] Adachi | Time-domain simulation of sound production in the brass instrument[END_REF]) This means in particular that the cost function lacks convexity. One typical solution to remedy this problem is to convexify the cost function using Tikhonov regularization, which amounts to adding quadratic terms with respect to some mask parameters. More precisely, we define the complete cost function C

M ref M Difference F ℓ (
C (p ref , M) = C • (p ref , p M ) + β Q Q 2 ℓ + β H H 2 (2.7)
The choice on the specific penalization has been made on Q ℓ and H because it proved to be • sufficient to have a well defined solution up to a sufficiently good precision (cf. table 2),

• necessary to remove very different solutions (cf. section 2.8).

It should be noted that this particular choice of penalization, instead of a more general form like (H -H 0 ) 2 for a reference H 0 which should be fixed for the whole optimization, implies that the optimization procedure will favor solutions with the smallest quality factor and lip opening. This was chosen for lack of a good candidate for H 0 .

This particular choice of penalization proved to give results close to those in the literature, except for H (cf. section 3). The method for fixing the values of β Q and β H is done so that a typical value of the penalization should be of the same magnitude as C • (p ref , p). As we expect C • (p ref , p) to be about 0.3 (cf. section 3), that Q ∼ = 7 and H ∼ = 10 -4 m (chosen among the known values in [START_REF] Velut | How Well Can Linear Stability Analysis Predict the Behaviour of an Outward-Striking Valve Brass Instrument Model?[END_REF]), we took

β Q = 5 × 10 -3 and β H = 3 × 10 7 m -2 .

Continuity, optimization algorithm

The algorithm chosen to find the minimum of the cost function is the dual annealing optimization (cf. [START_REF] Xiang | Generalized simulated annealing algorithm and its application to the Thomson model[END_REF]) which is a stochastic algorithm that requires neither the cost function to be regular, nor the minimum to be unique.

Indeed, the cost function in this article is not continuous: a small variation of the mask parameters can lead to completely dissimilar solutions. For example the trombone player can obtain the different notes by only varying its lip resonance frequency: the variation of playing frequency with lip resonance frequency clearly has jumps (cf. [START_REF] Mattéoli | Minimal blowing pressure allowing periodic oscillations in a model of bass brass instruments[END_REF] figure 4).

Another problem of the cost function is that it has many local minima. Although we do not have a mathematical or musical reason for this, it is clearly seen during the optimization process using dual annealing as it performs local searches before jumping to other locations (see table 3 where each line represents a local minimum).

Limitations of dual annealing

The dual annealing algorithm is known to give a solution for some very slow (logarithmically) decreasing temperatures (cf. [START_REF] Xiang | Generalized simulated annealing algorithm and its application to the Thomson model[END_REF]), but that kind of evolution of the temperature implies a very slow convergence. In practice, a faster decreasing temperature is used, but the convergence is not assured.

Moreover, as this algorithm is of stochastic nature, there is no simple criterion to stop it, and an arbitrary condition has to be chosen: the choice was made to bound by 1000 the number of calls to the cost function. With this choice, a typical run lasts about 1 hour on a desktop computer.

The probabilistic nature of this algorithm also means that two runs of the same algorithm, with the same initialization data (except for the seed of the random generator) give different solutions that can be quite different, as the global minimum might not have been reached in one run. It is therefore often necessary to launch the algorithm iteratively until a new run does not produce a solution with lower cost function. We say that a set of mask parameters satisfying this hypothesis is C -admissible. A random set of mask parameters is not C -admissible in general . Indeed, two sets M 1 and M 2 can give very similar waveforms (cf. table 1), so that C

Precision of the algorithm

• (p ref , p M1 ) = C • (p ref , p M2
), but have different quality factor or lip opening, and implying for example

C (p ref , M 1 ) < C (p ref , M 2 ). In that case, M 2 cannot be C -admissible.
Moreover, this definition highly depends on the choices made for the definition of C , be it α F or the choice of penalizations. A set of mask parameters may be C -minimal for one choice, but no longer for another one!

The assumption in this article is that for every "realistic" signal (i.e. coming from the recording of a trombone), there is only one mask parameter that is C -admissible. It is not at all clear that this is true, and this is even known to be false if the penalizations are not added (see section 2.5). Taking this assumption for granted, a set of mask parameters obtained by minimization of the cost function is automatically C -admissible.

To find a suitable C -admissible set of mask parameters and to stay close to an actual trombone signal, so that the robustness of the optimization procedure can be tested, we applied the optimization procedure to a reference signal on a recorded D4 (cf. section 2.2), cf. figure 1. The obtained values are given in the second column of table [START_REF] Avanzini | Modelling the Mechanical Response of the Reed-Mouthpiece-Lip System of a Clarinet. Part I. A One-Dimensional Distributed Model[END_REF], denoted by M ref . The set of mask parameters used for the initialization of the dual annealing algorithm is given in the third column, and the result of the optimization algorithm is in the last one. The search space is given in its caption.

M ref

Init The resulting waveforms for the two sets of mask parameters are indistinguishable, with a relative RMS error of only 3%, and an error on the frequencies of about 6 cents. The results of optimization are acceptable (much lower than the dispersion of the values found in the literature) and representative of the errors we found on other simulations (cf. section 3) except for the opening at rest H, but its value is so small that it is hard to give a physical interpretation (see section 3.2). The difference between M ref and M optim may be explained by the fact that we had to stop the algorithm at one point (cf. section 2.7) or because C is insufficiently convexified.

Continuation

During the analysis of the bend (cf. section 3.4), the continuation software auto-07p ( [START_REF] Doedel | Auto-07p, continuation and bifurcation software for ordinary differential equations[END_REF]) is used to follow the evolution of the playing frequency with different parameters. The continuation is first initialized for a reference set of mask parameters M 1 which is chosen for each musician to be the point which minimizes the cost function C • among all the optimized values of mask parameters, so as to be as close as possible to the actual recording.

The bifurcation diagram is built using the dependency on p m up to the recorded value of the reference signal, so that auto-07p is now precisely set to the signal p M1 .

Then a continuation curve along one of the mask parameters (either Q ℓ , F ℓ , H or µ), δ is computed, all other physical variables being fixed, and the playing frequency is drawn.

Results

The results obtained for the different musicians are presented in the following subsections, but first it is interesting to look at the mouth pressure as a function of the playing frequency, cf figure 2. Indeed, although all these data are directly recorded, and not optimized, we can clearly see differences between the two recording sessions of musician A, where the second session has a larger mouth pressure, which could translate into perceptible differences within the optimized data of a single musician.

Ab4

Errors on sustained notes

Both RMS and frequency errors obtained at the end of optimization are presented in figures 3 and 5. The RMS error can be quite large for some notes (up to 40%), which is not surprising as the model is one of the simplest and many physical details are neglected. As the optimization looks for a best fit among all mask parameters, this means the model should be complexified to take into account more of the physics of the instrument if precision on timbre and playing frequency have to be maintained, provided the dual annealing algorithm give results close enough to the global minimum. For reference, a typical waveform is shown in figure [START_REF] Bromage | Open Areas of Vibrating Lips in Trombone Playing[END_REF] where the reference signal is given in green, the reconstructed signal is in orange and the difference between them is in dotted red. The relative RMS error for this particular signal is 0.28. Although many properties are well approximated, the higher harmonics of the signal are clearly not in agreement with the experimental signal. This gives a typical value that can be expected for the RMS error. Concerning repeatability, estimations of the mask parameters of musician A are coherent and give almost the same results for both the RMS error and the frequency error (except for the playing frequency of the note Bb4). However, the errors differ largely between both players, player B mainly getting the lowest error. This may suggest that both musicians use different techniques, and that player B is closer to the simple model (2.1).

Note in particular the difference between errors for the note F3, where musician A has the largest RMS error of all (cf. figure 4), and musician B has one of the lowest.

Discussion on sustained notes

Lip resonance frequency

The lip resonance frequency as a function of the playing frequency is shown in figure 6 for all three musicians. As for any outward model, the frequency of the lips is lower than the playing frequency (cf. [START_REF] Campbell | Brass Instruments As We Know Them Today[END_REF]), which is clearly seen in this figure as the circles are below the line F play = F ℓ . Note that for a given frequency, there is little dispersion from player A (A1 or A2) to player B. Moreover the regression line

F ℓ = 0.9366F play -31.57 (3.1)
gives a good fit with R 2 = 0.996 and could be used as a first estimation of the playing frequency using only the lip frequency.

Quality factor

The estimated quality factors for all three players are displayed in figure 7. As in the case of lip frequency (cf. 3.2) results are very close for all three players, and also for all notes, being between 2 and 5. Compared to the literature, they are however smaller than the measured values of [START_REF] Cullen | Brass Instruments: Linear Stability Analysis and Experiments with an Artificial Mouth[END_REF] (between 9 and 10.5) but comparable to the estimation of [START_REF] Lopez | Physical Modeling of Buzzing Artificial Lips: The Effect of Acoustical Feedback[END_REF] (around 5), [START_REF] James | Experimental Mechanical and Fluid Mechanical Investigations of the Brass Instrument Lip-reed and the Human Vocal Folds[END_REF] (between 1.2 and 1.8), [START_REF] Richards | Investigation of the lip reed using computational modelling and experimental studies with an artificial mouth[END_REF] (around 3.7), [START_REF] Rodet | Physical models of trumpetlike instruments. detailed behavior and model improvements[END_REF] (around 2.88) and [START_REF] Adachi | Time-domain simulation of sound production in the brass instrument[END_REF] (between 0.5 and 3). Except for the first reference, this justifies the penalization on Q ℓ , which tends to favor the smallest possible values.

Ab4

The values obtained for the two recordings of player A are always very close, this may mean that it depends very little on the loudness.

Surface density

The optimized values of µ -1 are given in figure 8. Except for the 4 highest values, they are comparable to [START_REF] Rodet | Physical models of trumpetlike instruments. detailed behavior and model improvements[END_REF]. However, they are overestimated compared to other values found in the literature ( [START_REF] Elliott | Regeneration in brass wind instruments[END_REF][START_REF] Cullen | Brass Instruments: Linear Stability Analysis and Experiments with an Artificial Mouth[END_REF][START_REF] Lopez | Physical Modeling of Buzzing Artificial Lips: The Effect of Acoustical Feedback[END_REF][START_REF] Richards | Investigation of the lip reed using computational modelling and experimental studies with an artificial mouth[END_REF]) where they are between 0.03 and 0.2m 2 .kg -1 .

We can see that the data for A1 is systematically higher than that of A2, indicating a possible dependency on the mouth pressure and the loudness.

Opening at rest

The values of opening at rest obtained by optimization are given in figure 9. They seem very small compared to what was obtained by other authors, up to a factor 10: a typical value obtained by optimization is around 2.10 -5 m (see figure 9) whereas [START_REF] Cullen | Brass Instruments: Linear Stability Analysis and Experiments with an Artificial Mouth[END_REF] and [START_REF] Richards | Investigation of the lip reed using computational modelling and experimental studies with an artificial mouth[END_REF] have a typical value of 5.10 However, when shifting from opening at rest to mean opening (cf. figure 10) using formula (A.3), which takes also into account the mouth pressure, the lip frequency and the lip surface density, the results are comparable to those of [START_REF] Bromage | Open Areas of Vibrating Lips in Trombone Playing[END_REF], which are between 0.6mm and 2mm. It should be noted that in this case, the value of opening at rest is negligible in the formula (A.3) in appendix A.

Just as for µ -1 there appears to be a correlation between loudness and mean opening. It is not only expected, but actually obvious from the formula (A.3) where the mouth pressure appears.

Optimization for bent notes

The musicians were instructed to perform pitch bends on F4: without moving the slide, the player used embouchure adjustments to vary the pitch, first below its normal value, then above, then below, and then back to F4. For each recording of approximately 10 seconds, the signals are cut into chunks of 0.2 seconds, and the optimization procedure is applied independently on each chunk. The note F4 has been chosen because it is one of the most comfortable to bend for the musician. The RMS and frequency errors can be found in figures [START_REF] De | YIN, a fundamental frequency estimator for speech and music[END_REF] and [START_REF] Doedel | Auto-07p, continuation and bifurcation software for ordinary differential equations[END_REF]. The RMS error is around 0.25, except for a group of notes played by musician B with low frequencies, which may be related to a particular technique used by this musician. The frequency error is quite low except for the lowest notes, and is in agreement with the very low frequency error found for F4 in figure [START_REF] Campbell | Brass Instruments As We Know Them Today[END_REF]. This suggests that the model is not able to predict precisely what the musician is doing for the lowest frequencies of the bend. Indeed, there seems to be two regions for the errors, with a change at around 352Hz, as if the optimization process could not find parameters that fit the playing frequency below that value. In practice, musicians often use special techniques to bend to very low notes, such as using the vocal tract. This technique is clearly not taken into account in the model, and it seems the algorithm indicates its own limits.

Discussion on bent notes

During bending, the musician varies many parameters. This makes it quite difficult to see the influence of any of them. In the following diagrams, the evolution of playing frequency is shown with respect to the mask parameters. To put it into perspective, the theoretical evolution with respect to only the considered parameter (the other parameters being kept constant) is also computed using auto-07p (see section 2.9). The mask parameters used to initialize the continuation are those with smallest cost function among all the optimized values for this recording, to ensure that the model is as close to the measures as possible.

Quality factor Q ℓ

The results of the optimization for bent notes is presented in figure 13 for the quality factor. The values obtained are within the same range as in figure 7. One striking feature is the proximity from the measures of musician B, and the results of continuation obtained by auto-07p. It seems like the playing frequency is completely predicted by the evolution of the quality factor. However, the precision of the fit must be put into perspective with the rather large errors in the optimization (see figures 11 and 12). The fit is not so good with musician A, although the results of the continuation go in the right direction.

F4

Figure 13: Variation of Q ℓ as a function of F play . Crosses represent the measures for all three musicians, circles the initialization parameter for continuation (see section 2.9), and lines the continuation obtained with auto-07p. The dashed line represents the average frequency of the actual played F4 with no bend.

Lip resonance frequency F ℓ

The results of optimization for the lip resonance frequency are given in figure 14, and are quite difficult to interpret. Even more than in figure 12, there seem to be two regions, one before 352Hz, and one after.

Above 352Hz, the estimation of lip frequency does not give a clear tendency. Although we could expect the lip frequency to increase with playing frequency, just as in figure 6, this is not what appears in the figure. This suggests that the lip frequency is only a coarse tuner, and the quality factor is actually the fine tuner.

Lip surface density µ

The results of the optimization for bent notes is presented in figure 15 for the lip surface density. The values are compatible with those in figure 8, and the evolution of playing frequency relatively to µ is compatible with the theoretical one obtained by continuation.

Opening at rest H

The results of the optimization for bent notes is presented in figure 16 for the opening at rest. The values obtained are within the same range as in figure 9.

As explained in section 3.2, the values obtained for H are very small, and therefore not very well defined (see error term in table 2). The mean opening obtained from other optimized values and formula (A.3) is given in figure 17, and a clear tendency can be observed above 350Hz: the mean opening increases with the playing frequency for all musicians.

Below 350Hz the tendency is not so clear. Moreover, one must be careful with interpretation as there may be other phenomena involved than those directly modeled (cf. section 3.3). Crosses represent the measures for all three musicians, circles the initialization parameter for continuation (see section 2.9), and lines the continuation obtained with auto-07p. The dashed line represents the average frequency of the actual played F4 with no bend.

Optimization for a crescendo

The same procedure as in section 3.3 was used for the recordings of a crescendo on F4 for all three musicians.

The relative RMS error is presented in figure 18, and indicates that the higher the mouth pressure, the higher the RMS error. This shows that the simple model (2.1) is good at reproducing the timbre for low pressure, but not so much for higher pressures. This may be due to the nonlinear propagation along the length of the trombone. The error in frequency is presented in figure 19 and is compatible with that in figure 5. It proves that the model (2.1) is actually quite good at reproducing the frequency, whatever the dynamic of the playing. This indicates that the limits of the model are not so much on the frequency, but more on timbre. The waveform for two different dynamics are shown in figure 20, both for the measured signal, together with the reconstructed signal from optimized mask parameters. The difference in timbre is clearly seen for the forte recording. The details of the figures obtained through optimization during a crescendo can be found at Comparison of measured vs. simulated sounds with these parameters can be found at http://perso.univ-lemans.fr/~smauge/ mask.

Conclusion

In this article, a new method is proposed to estimate the mask parameters of a brass musician within a set of acceptable parameters (so called C -admissible). This approach is used on recordings of actual musicians during playing, and is able to deliver a coherent set of parameters (except maybe for the opening at rest), in that they are not too far from existing results in the literature, and their values evolve in a way that is compatible with theory during the playing. The values obtained prove that a simple model is already capable of reproducing a playing frequency close to that played by an actual musician, with a dynamic and waveform that are similar to the measured ones. This may prove useful for instrument making, although more research should be done to assess the robustness of the method, and investigate the variability of the mask parameters from player to player. Moreover, it can give new leads to a better understanding of intonation, such as the almost linear relation between playing and lip frequencies, or the role of quality factor as a fine tuner. Furthermore, the system seems to be able to detect when a particular technique is used, for example on the lowest part of bent notes where the vocal tract is used by the musician, as the difference between the techniques is clearly seen on the cost function.

  For a given set of mask parameters M ref , the goal is to find an algorithm that captures an approximation of M ref with only the knowledge of p ref and F ref . The algorithm we propose here is not able to do that without a prior assumption on M ref , namely (⋆) M ref must be a minimum of the cost function M → C (p ref , M)

Figure 1 :

 1 Figure 1: Block diagram explaining how to get a Cadmissible set of mask parameters from a recorded signal and use it to assess the precision of the algorithm.

Figure 2 :

 2 Figure 2: Measured mouth pressure P m averaged over one period as a function of the playing frequency for all three musicians and the six notes.

Figure 3 :

 3 Figure 3: RMS error of the signals obtained from optimized mask parameters, relative to the measured signal for 6 different notes and the three musicians.

Figure 4 :

 4 Figure 4: Typical waveform for the recorded signal (green) and for the signal obtained from optimized mask parameters (orange) as played by musician A2 on a Bb3. The difference between signals is drawn in dashed red. The relative RMS error is 0.28

Figure 5 :

 5 Figure 5: Frequency error in cents of the signals obtained from optimized mask parameters, relative to the recorded signal for 6 different notes and the three musicians

Figure 6 :

 6 Figure 6: Variation of F ℓ as a function of F play for all three musicians in circles, together with the diagonal F play = F ℓ in dashed black, and the regression line in dashed red.

Figure 7 :

 7 Figure 7: Variation of Q ℓ as a function of F play for all three musicians in circles.

Figure 8 :

 8 Figure 8: Variation of µ -1 as a function of F play for all three musicians.

Figure 9 :

 9 Figure 9: Variation of H as a function of F play for the three musicians.

Figure 10 :

 10 Figure 10: Mean opening of the lips as a function of F play for the three musicians.

F4Figure 11 :

 11 Figure 11: RMS error for bent note on F4 for all three musicians. The dashed line represents the average frequency of the actual played F4 with no bend.

F4Figure 12 :

 12 Figure 12: Frequency error for bent note on F4 for all three musicians. The dashed line represents the average frequency of the actual played F4 with no bend.

Figure 14 :

 14 Figure14: Variation of F ℓ as a function of F play . Crosses represent the measures for all three musicians, circles the initialization parameter for continuation (see section 2.9), and lines the continuation obtained with auto-07p. The dashed line represents the average frequency of the actual played F4 with no bend.

F4Figure 15 :

 15 Figure15: Variation of µ -1 as a function of F play . Crosses represent the measures for all three musicians, circles the initialization parameter for continuation (see section 2.9), and lines the continuation obtained with auto-07p. The dashed line represents the average frequency of the actual played F4 with no bend.

F4Figure 16 :

 16 Figure16: Variation of H as a function of F play . Crosses represent the measures all three musicians, circles the initialization parameter for continuation (see section 2.9), and lines the continuation obtained with auto-07p. The dashed line represents the average frequency of the actual played F4 with no bend.

Figure 17 :

 17 Figure 17: Variation of H mean as a function of F play .Crosses represent the measures for all three musicians, circles the initialization parameter for continuation (see section 2.9), and lines the continuation obtained with auto-07p. The dashed line represents the average frequency of the actual played F4 with no bend.

Figure 18 :

 18 Figure 18: Variation of RMS error as a function of P m for all three players.

Figure 19 :

 19 Figure 19: Variation of frequency error as a function of P m for all three players.

Figure 20 :

 20 Figure 20: Comparison of wave functions for measured signal vs. optimized signal for musician B. Top for piano (relative RMS error : 20%, playing frequency error 1.5 cents)), and bottom for forte RMS error : 34%, playing frequency error 0.7 cents))

F

Table 1 :

 1

	Hz)	177.70	184.05	60 cents
	Q ℓ	3.82	4.79	20%
	µ (kg.m -2 )	1.28	1.42	10.9%
	H (m)	1.2 × 10 -5 1.9 × 10 -4	1483%

Values of two different sets of mask parameters giving almost identical signals: relative RMS difference is 0.9% and difference in frequencies is 1.41 cents. Mouth pressure is equal for both simulations and fixed at 2500Pa.

Table 2 :

 2

	.	M optim (error)

Data for optimization: Mouth pressure p m = 1656P a and width of the lips w = 12.10 -3 m. Search space is F ℓ ∈ [150, 200], Q ℓ ∈ [0.1, 6], µ ∈ [0.1, 3], H ∈ [10 -5 , 10 -3 ]. Error on playing frequency : 6cents

Table 3 :

 3 ℓ (Hz) Optimization steps. Each line represents the result of a gradient descent performed by the dual annealing algorithm Note F ℓ (Hz) Q ℓ µ(kg.m -2 ) H(mm) P m (P a) F play (Hz)

		Q ℓ		µ(kg.m -2 ) H(mm)	Cost	Frequency (Hz)
	184.38090 4.27669	1.27669	0.09	0.11377	233.18
	184.33001 4.54898	1.04458	0.08	0.10951	233.18
	180.35897 4.34470	1.14597	0.11	0.10177	233.04
	178.16236 4.37235	1.17461	0.15	0.09876	232.95
	176.77990 4.35483	1.02964	0.09	0.09696	232.98
	176.77899 4.33623	1.03454	0.11	0.09614	233.01
	175.21228 4.10074	1.10552	0.11	0.09348	233.04
	175.21124 4.10265	1.09939	0.13	0.09318	233.08
	175.21124 4.10265	1.09939	0.03	0.09000	232.93
	175.21816 4.08641	1.05548	0.01	0.08922	232.97
	174.56073 4.13190	1.07071	0.03	0.08912	232.91
	174.56073 4.13190	1.07071	0.03	0.08909	232.91
	174.56073 4.15607	1.07071	0.03	0.08718	232.71
	174.56169 4.15607	1.07071	0.03	0.08613	232.47
	174.56169 4.15610	1.07071	0.03	0.08608	232.43
	174.56160 4.15610	1.07071	0.03	0.08608	232.40
	174.56160 4.15610	1.07015	0.03	0.08600	232.56
	174.56098 4.15610	1.07015	0.03	0.08593	232.52
	174.55610 4.15610	1.07015	0.03	0.08587	232.47
	174.55610 4.15539	1.07015	0.03	0.08548	232.49
	Bb2	78.78	3.18	1.33	0.32	896	116.12
	Bb3	167.71 3.54	0.57	0.0114	1433	234.26
	Bb4	392.90 3.42	0.62	0.452	5134	469.69
	D4	232.11 4.56	1.32	0.38	2488	295.62
	F3	123.77 2.40	1.22	0.0139	1350	175.45
	F4	277.53 3.79	1.68	0.0132	3716	351.19

Table 4 :

 4 Values obtained by optimization for different notes played by musician B (mezzo forte) for a width w = 0.012m. Comparison of measured vs. simulated sounds with these parameters can be found at http://perso.univ-lemans. fr/~smauge/mask/#sounds

	s n	C n
	-12.24 + 237.80i	237.83 + 12.23i
	-16.01 + 700.45i	210.28 + 4.81i
	-21.08 + 1063.24i	222.43 + 4.41i
	-24.09 + 1438.12i	210.97 + 3.53i
	-26.14 + 1838.74i	237.99 + 3.38i
	-30.00 + 2176.32i	268.08 + 3.70i
	-32.85 + 2532.56i	159.47 + 2.07i
	-35.80 + 2920.16i	150.19 + 1.84i
	-40.08 + 3307.88i	175.93 + 2.13i
	-45.82 + 3707.10i	99.70 + 1.23i
	-276.05 + 4886.60i 1036.27 + 58.54i

Table 5 :

 5 Values of the coefficients in equation for the impedance decomposition for the Bb trombone in first position 2.2

Acknowledgment

The authors would like to thank Christophe Vergez for many discussions and very helpful comments on the manuscript.

In memoriam

This article is dedicated to the memory of Joël Gilbert , who was at the origin of this project and a constant source of ideas.

A Mean opening

Suppose that p and h are of period T , write for the mean value of any T -periodic function. Applying it on equation

gives

As h is T -periodic, h ′ = h ′′ = 0, and as p = 0 we get

B Tables