© A. Ernoult et al., Published by EDP Sciences, 2021

1 Introduction

A wind instrument can be divided into two parts: its exciter which creates an acoustic oscillation under the musician’s action, and its resonator in which acoustic waves propagate. The resonator is essentially a tube, characterized by the variation of its internal cross-sectional area along its main axis and by the presence of side holes with specific locations and dimensions. In this study, the bore profile encompasses all these geometric parameters. In traditional instrument organology, most woodwind instruments have side holes, unlike most brass instruments. By extension, in this study, woodwind-like instruments refer to instruments with side holes.

The bore reconstruction of such instruments aims at obtaining their bore geometry from acoustic measurements. This non-destructive characterization can be useful for exploring historical instruments, for documenting museum collections or for building reproductions of them. It can also be used for detecting defects on built instruments. This technique can finally be seen as a first step towards automating the design of new musical instruments by taking a somewhat similar approach but using targeted acoustic specifications instead of acoustic measurements.

In the literature, bore reconstructions from measured data have only been performed on instruments without side holes (brass or section of woodwind instruments with no open side holes). The first approach was using pulse reflectometry [1, 2] in which the variation of the cross sectional area along the main pipe axis is computed directly from the input impulse response (in the time domain). The reconstructed bore takes the form of a stepped tube, the axial resolution being determined by the highest measured frequency. This method has been used to reconstruct trumpets [1, 3] or bassoon crooks [4, 5]. This approach cannot be directly applied to instruments with side holes where effectively the waves propagate in several tubes in parallel. However, reflectometry has been used by Chilekwa [6] to locate hole and by Mamou-Mani et al. [7] to investigate the difference between holes on different oboes. In this last case the reconstructed geometry obtained for the upstream main bore section with closed holes was not directly comparable to the actual geometry due to the propagation within the closed chimneys.

As a second approach, optimization algorithms can be used to find the geometry that minimizes the deviation between simulations and acoustic measurements. Kausel [8] minimizes the residuals between target and simulated impedances of a brass instrument by using a zeroth order algorithm (Rosenbrock). This type of local algorithm, in which no gradient is computed, is robust but not efficient in terms of convergence rate and computational cost. A “virtual trumpet” with up to 100 design variables is reconstructed from simulated targets by imposing the mouthpiece geometry, and an elementary stepped tube is reconstructed from measurements. In a slightly different approach, Chilekwa [6] uses the Rosenbrock algorithm to localize and characterize defects in a cylindrical duct, which can be considered as the reconstruction of a really basic woodwind instrument. Schmutzhard et al. [9] proposed a particular expression in the time domain for thermo-viscous losses. He used optimization to obtain the shape of the bore of such a virtual instrument which gives an acoustic response equivalent to that measured on a real instrument. He adjusted up to 20 cones by minimizing the L²-norm of the difference between the impedances with a gradient-based algorithm (Levenberg–Marquardt) using finite difference to compute the gradient.

As mentioned previously, the automated design of instruments can be seen as a similar issue. For this purpose, Noreland [10] used a gradient based algorithm with analytic gradient to design a brass bore profile composed of up to 400 truncated cones, by optimizing the frequency and magnitude of resonances. This type of first order algorithm is more efficient but requires the gradient of the cost function to be estimated. Braden et al. [11] tuned some properties of the impedance of a trombone by using the Rosenbrock algorithm to minimize three cost functions relating to the total impedance, the frequency and the magnitude of the resonances. He parameterized the geometry of the brass bell in terms of a few Bessel horns to reduce the number of design variables. Colinot et al. [12] designed a coaxial saxophone (without side hole) equivalent to an existing conical instrument by minimizing the p-norm of the impedance difference.

Focusing on woodwind instruments, Debut et al. [13] presented one of the first attempts to design instrument with side holes using an optimization process. Based on length correction considerations, he adjusted the location of the register hole of a clarinet to improve the tuning of the second register. He also tried to optimize some characteristics of the main bore in the same purpose. Lefebvre [14] designed simple flute-, clarinet-, or saxophone-like instruments with up to seven tone holes. He used the L-BFGS-B algorithm with gradients estimated by finite-differences, to minimize the deviation between the first resonance frequency and the expected playing frequency for several fingerings. The design variables were the location and the radius of the holes and the main bore length. With a similar approach, Noreland [15] followed by Guilloteau et al. [16] designed a logical clarinet (one hole per note without any fork fingering) by optimizing the resonance frequencies of the two first registers, then the ratio of their magnitudes. They used a similar gradient-based algorithm to adjust the dimensions of the holes (location, radius, chimney height), the main bore geometry being fixed. In all these studies, the gradient is computed by finite-difference, due to the apparent complexity of its explicit computation in the presence of holes. This approximation can lead to numerical issues as illustrated by Ernoult [17] who designed a more complex instrument in which the main bore was also optimized. He showed that the traditional resonance definition leads to a discontinuous cost-function which cannot be easily optimized. A new and more robust definition is proposed to tune the resonance frequency and magnitude. In a slightly different approach, Tournemenne et al. [18] proposed tuning some geometric parameters of a brass instrument (up to 10) by minimizing a cost function based on sound characteristics averaged on simulated sounds, derived from a set of musician parameters. He used a surrogate-assisted derivative-free optimization strategy, particularly suited to this non differentiable problem.

The aim of the present study is to find a robust and efficient methodology to reconstruct instruments with side holes from acoustic measurements. Following the literature, the reconstruction is performed by minimizing the deviation between simulations and measurements of an observable derived from the input impedance for several fingerings. To solve this inverse problem, computing the gradient of the cost function without using finite differences leads to a faster and more accurate estimation. This reconstruction process being a first step towards instrument design, a versatile formulation in which the cost function can be easily changed is developed. The full waveform inversion (FWI) appears perfectly suited to this purpose. This technique, which has its origins in the seismology community [19], is based on the observation that the gradient of the acoustic fields (here flow and pressure) in each point of the medium with respect to the design variables satisfies the same equations as the forward problem of modeling the wave propagation in a known propagation medium, just with different source terms. This enables the gradient estimation of the inverse problem with a negligible additional computation cost and a better accuracy compared to the finite difference method used in most previous studies.

The wave propagation model and the finite element method used to compute the impedance (forward problem) through the inversion of a linear system are presented in Section 2. In Section 3, the general bore reconstruction problem (inverse problem) is presented. The FWI is applied to woodwind instruments, allowing an explicit characterization of the gradient of the cost function. Considerations regarding the observable are carefully laid out for the case of an elementary instrument in Section 4, resulting in a reconstruction strategy applicable to measured impedance data. The performance of this technique is then discussed, and a conclusion is given in Section 5. All the tools presented here are implemented and available in the open source software OpenWInD [20] under GPLv3.

2 Forward problem

2.1 Instrument description

The reconstruction methodology (theory and strategy) is presented with reference to a real resonator with a simple geometry for which it is possible to measure manually the design parameters (Fig. 1). The main bore (index mb in the rest of the document) is a cylindrical pipe seen, in the reconstruction problem, as a special example of a conical duct, described by three design parameters: its total length L_mb = 287.5 mm and the radii at the left and right ends, both equal to R_mb (0) = R_mb (L_mb) = 2 mm. It is drilled with four side holes (h1, h2, h3, h4) forming cylindrical chimneys, each characterized by three design parameters: their location on the main pipe X_hi, their chimney height H_hi and their radius R_hi. The values obtained by direct measurement with caliper or graduated scale are indicated in Table 1 with the corresponding estimated uncertainties (the chimney heights being measured indirectly, their relative uncertainties are higher than other parameters).

Figure 1 Sketch of the cylindrical tube with four side holes to be recovered. It is described by 15 geometric parameters (three for the main bore profile and three for each of the four holes).

The instrument geometry is therefore described by N_ξ = 15 design parameters grouped together in the vector ξ. The forward problem presented in the current section, allows the computation of the acoustic response of the instrument in the frequency domain, including its input impedance Z, from the knowledge of these design parameters.

To model the acoustic propagation, the instrument is described by a network of pipes and connectors which defines the pipe’s boundary conditions (Fig. 2). The studied instrument is composed of 9 pipes and 10 connectors: 5 pipes for the main bore and one for each of the four hole chimneys. In these pipes, the wave propagation is modeled in one dimension by telegrapher’s equations including thermo-viscous effects, detailed in Section 2.2. At the hole locations, the three pipes (upstream, downstream and chimney) are connected by a T-joint junction including acoustic masses, detailed in Section 2.3. The open end of each pipes is modeled by a radiation condition, developed in Section 2.4. As an input condition, at each frequency, a unitary flow is imposed at the entrance of the instrument.

Figure 2 The network representation of the studied resonator is composed of 9 pipes and 10 connection conditions.

2.2 Acoustic propagation along the pipe

Each pipe n ∈ [1, N] is characterized by the change in radius R_n (x, ξ) along its main axis x and its length L_n (ξ), both depending on the design parameters grouped in the vector ξ. The wave propagation is modeled by telegrapher’s equations, following the methodology presented by Tournemenne and Chabassier [21]. The variations of the acoustic pressure p_n and of the acoustic flow u_n along the nth pipe of length L_n, at the angular frequency ω, follow the equations ([22], Chap. 5.5, p. 239):

$$\{\begin{array}{ll}{Z}_v(x,\omega )u_n+\frac{\mathrm{d}{p}_n}{\mathrm{d}x}=0& \\ Y_t(x,\omega )p_n+\frac{\mathrm{d}{u}_n}{\mathrm{d}x}=0& \end{array} \hspace{1em}\forall x\in \left[0,L_n\right],$$(1)

with thermo-viscous effects included through:

$$\{\begin{array}{c}{Z}_v\left(x,\omega \right)=\mathrm{jw}\frac{\rho }{S_n\left(x\right)}{\left(1-\mathcal{J}\left(k_v{R}_n\left(x\right)\right)\right)}^{-1}\\ Y_t\left(x,\omega \right)=\mathrm{jw}\frac{S_n\left(x\right)}{\rho c^2}\left(1+(\gamma -1)\mathcal{J}\left(k_t{R}_n\left(x\right)\right)\right)\end{array}$$(2)

$$k_v=\sqrt{\mathrm{j\omega }\frac{\rho }{\mu }}, \hspace{1em}{k}_t=\sqrt{\mathrm{j\omega }\frac{\rho C_p}{\kappa }},$$

where S_n(x) = πR_n(x)² is the cross sectional area at the position x and $\mathcal{J}:\mathbb{C}\mapsto \mathbb{C}$ is a complex function modeling the viscothermal terms ([22], Chap. 5.5, p. 239)

$$\mathcal{J}\left(z\right)=\frac{2}{z}\frac{J_1\left(z\right)}{J_0\left(z\right)}\hspace{1em}\forall z\in \mathbb{C}$$(3)

with J₀ (z) and J₁ (z) Bessel functions of the first kind. Table 2 describes the physical constants, which are assumed uniform on the entire instrument in this study.

In equation (1), each individual pipe length L_n is taken into account indirectly through the definition domain of x. To explicitly write these physical parameters in the equations, and provide the possibility modifying them without changing the definition domain, the scaling $x=L_n\tilde{x}$ is used:

$$\{\begin{array}{ll}{L}_n{Z}_v\left(\tilde{x},\omega \right)u_n+\frac{\mathrm{d}{p}_n}{\mathrm{d}\tilde{x}}=0& \\ L_n{Y}_t\left(\tilde{x},\omega \right)p_n+\frac{\mathrm{d}{u}_n}{\mathrm{d}\tilde{x}}=0& \end{array}\hspace{1em}\forall \tilde{x}\in \left[\mathrm{0,1}\right].$$(4)

Now, all the physical parameters are coefficients of these equations. With this substitution, the optimization process only modifies the coefficients of the equations with the definition domain and topology remaining unchanged. The boundary conditions at $\tilde{x}=0$ and $\tilde{x}=1$, necessary to complete this system, depend on the type of connection at each extremity of the pipe.

2.3 T-junction

At the location of a side hole, the pressure and acoustic outward flow at the boundary of the chimney pipe (p₃, −u₃), at the upstream (p₁, u₁) and downstream (p₂, −u₂) main bore parts (Fig. 3), are related through the linear system ([22], Chap. 7.7 p. 364–370)

$$\begin{array}{c}\begin{array}{cc}\mathrm{j\omega }\left(\begin{array}{cc}{m}_{11}& m_{12}\\ m_{12}& m_{22}\end{array}\right)& \left(\begin{array}{c}{u}_1\\ -u_2\end{array}\right)\end{array}=\left(\begin{array}{ccc}1& 0& -1\\ 0& 1& -1\end{array}\right)\left(\begin{array}{l}{p}_1\\ p_2\\ p_3\end{array}\right)\\ u_1-u_2-u_3=0\end{array}$$

where m₁₁, m₁₂ and m₂₂ are acoustic masses linked to the junction geometry. In the case of a symmetric junction, m₁₁ = m₂₂, which is assumed within this study. They are generally expressed through the masses m_s and m_a: m₁₁ = m₂₂ = m_s + (m_a/4) and m₁₂ = m_s − (m_a/4), which appear during diagonalization of the mass matrix. This leads to the introduction of two auxiliary variables γ₁₂₃, ζ₁₂₃, such that:

$$\{\begin{array}{c}\left(\begin{array}{c}{u}_1\\ -u_2\\ -u_3\end{array}\right)={-T}^T\left(\begin{array}{c}{\gamma }_{123}\\ {\zeta }_{123}\end{array}\right)\\ \mathrm{j\omega }M\left(\begin{array}{c}{\gamma }_{123}\\ {\zeta }_{123}\end{array}\right)+T\left(\begin{array}{c}{p}_1\\ p_2\\ p_3\end{array}\right)=0\end{array}$$(6)

with

$$M=\left(\begin{array}{ll}2m_s& 0\\ 0& m_a/2\end{array}\right), \hspace{1em}\mathrm{and}\hspace{1em} T=\frac{1}{\sqrt{2}}\left(\begin{array}{lll}1& -1& 0\\ -1& -1& 2\end{array}\right).$$(7)

Figure 3 Sketch of a T-junction with the adopted convention for the flow orientation.

For circular pipes, the masses are linked to the geometric parameters by polynomial expression fitted on data obtained via 3D simulation [23, 24]. Here, the following expression for m_s is assumed [23]:

$$m_s=\frac{\rho }{\pi R_h}\left(0.82-0.193\delta -1.09{\delta }^2+1.27{\delta }^3-0.71{\delta }^4\right)+m_{\mathrm{mv}},$$(8)

with δ = R_h/R_mb and $R_h,R_{\mathrm{mb}}$ respectively the inner radii of the chimney pipe and of the main bore pipe at the hole location. In the chosen representation, the matching volume mass [25]

$$m_{\mathrm{mv}}=\frac{\rho }{\pi R_h}\frac{\delta }{8}\left(1+0.207{\delta }^3\right),$$(9)

is included in the symmetric mass and not added to the effective length of the chimney pipe as usual [23, 24]. Indeed, despite the fact that both expressions are similar at the limit $\frac{\omega m_{\mathrm{mv}}c{S}_{\mathrm{hole}}}{\rho }\ll 1$ and that the wave propagation in this volume is debatable, this representation incorporates all the coefficients depending on both the main bore and the side hole geometry into the junction. This is consistent with the network representation (Fig. 2). This representation also prompts us to use the long chimney approximation for the acoustic mass m_a [23]

$$m_a=\frac{\rho R_h}{\pi R_{\mathrm{mb}}^2}\left(-0.37+0.087\delta \right){\delta }^2,$$(10)

in which it is independent of the open/close state at the other extremity of the chimney pipe. This approximation is valid when the chimney height is longer than the chimney radius ($H_h>R_h$), however the induced error stays small for the studied geometry (below 2% of deviation from the complete formula for the 4 holes). This approximate model suffers from the major drawback of being non-passive in some practical cases, and could be improved using asymptotic approaches as [26, 27].

2.4 Radiation condition

Pipe radiation is modeled by a radiation admittance Y_r,n at each opening (the bell and the 4 holes). For an opening of radius R, it follows that:

$$u_n\left(1\right)=Y_{r,n}(\omega ,R)p_n\left(1\right).$$(11)

The radiation admitances are computed with the “non-causal” formula given by Silva et al. [28]:

$$Y_{r,n}\left(\omega ,R\right)=\frac{\pi R^2}{\rho c}\frac{1}{j\tan \left(\mathrm{kR}\mathcal{L}+0.5j\log \left(\left|\mathcal{R}\right|\right)\right)},$$(12)

with k = ω/c the wavenumber and

$$\mathcal{L}=\eta \frac{1+b_1(\mathrm{kR}{)}^2}{1+b_2(\mathrm{kR}{)}^2+b_3(\mathrm{kR}{)}^4+b_4(\mathrm{kR}{)}^6},$$(13)

$$\left|\mathcal{R}\right|=\frac{1+a_1(\mathrm{kR}{)}^2}{1+(\beta +a_1)(\mathrm{kR}{)}^2+a_2(\mathrm{kR}{)}^4+a_3(\mathrm{kR}{)}^6}$$(14)

with η, β, a₁, a₂, a₃, b₁, …, b₄ real coefficients depending on the radiation condition [28]. The choice of this condition can impact the result of the reconstruction. It is assumed here that the holes radiate with an infinite flanged condition, and the main bore with an unflanged condition. This choice is discussed in more detail in Section 4.4. For closed holes, the flow is set to zero by assuming Y_r,n = 0.

2.5 Numerical resolution

To solve the overall problem, the discretization of equations (4), (6), (11) is carried out using the high order spectral finite elements method in 1D and setting a unitary flow at the pipe entry (input). Following Tournemenne and Chabassier [21], the pressure and flow inside each of the individual pipes within the instrument are sought as piecewise polynomials of variable order on small intervals called elements. The interpolation points and the quadrature points are chosen as the Gauss Lobatto points on each element, which include the end points. The junction and radiation equations relate the acoustical variables as intended. The number of discrete unknowns increases with the polynomial order and with the number of elements. The method converges exponentially with the polynomial order, achieving arbitrarily accurate approximations for low computational costs.

Practically, a linear system must be solved at each frequency:

$$\mathit{A}\left(\omega ,\mathbf{\xi }\right)\mathbf{U}=\mathbf{E},$$(15)

where $\mathbf{U}{=}^t(u_{h,1},p_{h,1},\dots , u_{h,N},p_{h,N},{\gamma }_{123},{\zeta }_{123},\dots ,$ ${\gamma }_{(N-2)(N-1)N},{\zeta }_{(N-2)(N-1)N})\in {\mathbb{C}}^{N_h}$ is the vector of the pressure and the flow at each degree of freedom of the mesh, along with the junctions’ auxiliary variables. The acoustic fields (flow and pressure) at any point of the instrument can be recovered from the vector of unknowns U. In the vector E, the degree of freedom corresponding to the entrance of the instrument is set as a unitary flow. The matrix A(ω,ξ) has a sparse block structure, related to the organization of the unknowns inside U. Equation (16) illustrates this structure for a pipe with a single side hole.

$$\mathit{A}=\left(\begin{array}{llllllll}{A}_1& 0& 0& 0& 0& & & \\ & & & -\frac{1}{\sqrt{2}}{C}_1& \frac{1}{\sqrt{2}}{C}_1& & & \\ 0& A_2& 0& 0& 0& & & \\ & & & \frac{1}{\sqrt{2}}{E}_2& \frac{1}{\sqrt{2}}{E}_2& & & \\ 0& 0& A_3& 0& 0& & & \\ & & & 0& -\sqrt{2}{E}_3& & & \\ 0& \frac{1}{\sqrt{2}}{C}_1^{*}& 0& -\frac{1}{\sqrt{2}}{E}_2^{*}& 0& 0& \mathrm{j\omega }2m_s& 0\\ 0& -\frac{1}{\sqrt{2}}{C}_1^{*}& 0& -\frac{1}{\sqrt{2}}{E}_2^{*}& 0& \sqrt{2}{E}_3^{*}& 0& \frac{\mathrm{j\omega }{m}_a}{2}\end{array}\right)$$(16)

where the “asterisk” indicate the trans-conjugate of the matrices and with:

$$A_1=\left(\begin{array}{ll}\mathrm{j\omega }{M}_1^u& -B_1\\ B_1^T& \mathrm{j\omega }{M}_1^p\end{array}\right),$$

$$A_n=\left(\begin{array}{ll}\mathrm{j\omega }{M}_n^u& -B_n\\ B_n^T& \mathrm{j\omega }{M}_n^p+Y_{r,n}\left(\omega \right)C_n\end{array}\right)\hspace{1em} \mathrm{for} n\in \left[2, 3\right].$$

The matrices $M_n^u,M_n^p$ and B_n follow from the finite element discretization of equation (4), the matrices C_n and E_n relate respectively to the downstream and upstream end of the pipes. The expressions of these matrices are detailed in [21]. They are sparse, and depend on the values of the design parameters ξ. For any set of design parameters ξ and angular frequency ω, equation (15) defines the solution U, so in a sense,

$$\mathbf{U}\equiv \mathbf{U}(\omega ,\mathit{\xi })={\left(\mathit{A}(\omega ,\mathit{\xi })\right)}^{-1}\mathbf{E}.$$(16)

In this study, only observables based on the normalized input impedance are analyzed. The acoustic flow being set to 1 at the entrance (see Sect. 2.1), the input impedance corresponds to the pressure at the entrance of the instrument. The normalized input impedance $Z \in \mathbb{C}$ is therefore obtained from the simulated acoustic field by using the restriction matrix $\mathcal{R}\in {\mathbb{R}}^{1\times N_h}$:

$$Z(\omega ,\mathit{\xi })\equiv Z\left(\mathit{U}\right(\omega ,\mathit{\xi }\left)\right)=\mathcal{R}\mathit{U}(\omega ,\mathit{\xi }).$$(18)

Due to our choice of finite elements, the restriction matrix has only one non-zero value. Here, with i₀ the index of the degree of freedom corresponding to the pressure at the entrance of the instrument:

$${\mathcal{R}}_i=\{\begin{array}{ll}\frac{S_1\left(0\right)}{\rho c}& \mathrm{if}\hspace{1em}i=i_0\\ 0& \mathrm{otherwise}\end{array}$$(19)

with S₀(0) the cross sectional area at the entrance of the instrument. The input impedance is computed for each of the N_ω frequencies and for each of the N_fing fingerings (N_ω × N_fing matrix inversions) using the python 3 toolbox OpenWInD [20].

2.6 Measurement-simulation agreement

For the bore reconstruction, the previously described model is used to estimate the geometry of the resonator from a measured impedance. To ensure the validity of the approach, simulations and measurement must be verified to match on a known geometry in the studied frequency range.

In this study, the impedances are measured with a home-made sensor based on the two microphones, three calibrations method proposed by Gibiat and Laloë [29]. The measurement apparatus comprises a 2 mm inner radius and 450 mm long cylinder equipped with 5 microphones (Audix TM1) located at [6, 18, 51, 174.8, 400] mm from the entrance, enabling the impedance to be measured between 50 Hz and 12 kHz. The data analysis approach proposed by Dickens et al. [30] is used, which allows measurements to be made with only two calibration steps (blocked: infinite impedance, 100 m tube: unitary impedance). In order to compare the measurements acquired at the temperature T_meas to the simulation computed at 20 °C, the change of the sound velocity is corrected by multiplying the frequency axis by $\sqrt{(20+T_0)/(T_{\mathrm{meas}}+T_0)}$. The effect of the temperature variation on the thermo-viscous losses is not corrected (<2% of deviation between the impedances computed at 20 °C and 25 °C after frequency correction).

The agreement between measurements and simulations is validated on a cylinder with an inner radius of 2 mm and a length of 436 mm, for which the theoretical impedance is well known ([22], Chap. 5.5, p. 239–248). The relative error ($\left|\right|Z_{\mathrm{theo}}-Z_{\mathrm{meas}}|{|}_{\mathcal{l}{l}^2}/|\left|Z_{\mathrm{meas}}\right|{|}_{\mathcal{l}{l}^2}$) in the frequency range [100 Hz, 4 kHz] is about 5%. This corresponds to a deviation within ±2 cents for the frequency of the peaks and within ±0.5 dB for their magnitude. This is consistent with the deviation obtained with other impedance sensors [31].

This deviation mainly comes from measurement uncertainties. The electric and acoustic noise in the acquired data contribute to it, but can be smoothed by averaging. The acoustic loads (blocked and 100 m tube) used for the calibration steps are not perfect due to visco-thermal effects at the wall and bends for the long tube. This results in a small bias in the measured impedance. The model used here (Sect. 2.2) assumes only one propagation mode corresponding to kR ≪ 1. The same assumption is made for the Zwikker and Kosten model (Eq. (2)). For the presented geometry, this corresponds to f ≪ 30 kHz. The model is therefore theoretically valid over the entire frequency range studied here. However the radiation condition used (infinite flanged or unflanged pipe) does not perfectly correspond to the experimental situation. This approximation also contributes to the deviation between measurement and simulation.

3 Inverse problem

The bore reconstruction is performed by solving an optimization problem. The objective is to find the values of the design parameters which minimize the deviation between the simulated impedances (Sect. 2.5) and the measured ones (Sect. 2.6). This optimization problem is formulated here for the elementary instrument presented in Section 2.1 and can be easily generalized to any other configuration.

3.1 Initial geometry

The method presented in this study is not able to modify the topology of the resonator which must be known and fixed at the initial state: here a conical main pipe with four side holes (Sect. 2.1). Only the inner radius at the entrance of the pipe is assumed to be known (R_mb 0 = 2 mm). The aim of the inverse problem is to obtain the value of the 14 other geometric parameters from the measured input impedances: the main bore length L_mb, the main bore right end radius R_mb (L_mb), the locations of the 4 holes X_hi, their chimney heights H_hi and their radii R_hi.

The initial values given to these design parameters are indicated in Table 3 and can be compared to the values of Table 1 which were measured directly using calipers.

3.2 Constraints and parameterization

The values of the geometrical parameters of woodwind instruments have strong physical limitations. They all must be strictly positive, the holes must be located on the main bore and their radii must be strictly smaller than the main bore radius at their location. The resulting bounds are summarized in Table 3 and are ensured by defining a bound-constrained minimization problem. These algorithms do not necessarily guarantee that the bounds are respected for each evaluation of the cost function. A margin of 10⁻⁷ is therefore included in the bounds but omitted for clarity in what follows. The positive nature of R_mb (L_mb), L_mb and H _hi is ensured by fixing their bounds to ξ_min = 0 and ξ_max = ∞. To guarantee the inequalities between the radii of the holes and the main bore radius and between the locations of the holes and the main bore length, the design parameters modified by the optimization algorithm are chosen as auxiliary parameters ξ_p and ξ_q defined as:

$$X_{\mathrm{hi}}={\xi }_p{L}_{\mathrm{mb}},$$(20)

$$R_{\mathrm{hi}}={\xi }_q{R}_{\mathrm{mb}}\left(X_{\mathrm{hi}}\right)$$(21)

and bounded by 0 and 1.

3.3 General optimization problem

The reconstruction of the instrument is carried out by solving a least-square optimization problem. The aim of the process is to minimize the difference between the simulated and the measured observable (extracted from the acoustic fields) by adjusting the values of the N_ξ design parameters. The simulated impedance is compared to the measured one through an observable:

$$\begin{array}{llll}\varphi :& \mathbb{C}& \to & \mathbb{C}\\ & Z& \to & \varphi \left(Z\right)\end{array}$$(22)

which can be the identity, the modulus, or any other non linear function. The choice of this observable is crucial to observe a correct convergence of the algorithm. This choice is discussed in Section 4.1 for the studied instrument. The only limitation on the choice of ϕ is its differentiability with respect to Z in the complex sense (see Appendix A). This ${\mathcal{C}}^1$ property is necessary for the use of gradient-based optimization algorithms (Sect. 3.4).

The problem is treated frequency by frequency and fingering by fingering. For a given fingering, and for each angular frequency ω, the real residual r(ω, ξ) is defined as

$$r\left(\omega ,\mathbf{\xi }\right)\equiv r\left(\varphi \left(Z\left(\omega ,\mathbf{\xi }\right)\right)-\varphi \left(Z_{\mathrm{meas}}\left(\omega \right)\right)\right),$$(23)

where

$$r\left(\zeta \right)=\left(\begin{array}{l}\mathfrak{R}e\left(\zeta \right)\\ \mathfrak{I}m\left(\zeta \right)\end{array}\right)\hspace{1em}\forall \zeta \in \mathbb{C}.$$

All the frequencies for all the fingerings are concatenated in the vector r(ω, ξ) of length 2 × N_ω × N_fing. For clarity, the ω dependency is omitted in the rest of the document. The choice of using a real residual is motivated by the gradient computation detailed in Section 3.4. The cost function is defined as:

$$F\left(\mathit{\xi }\right)=f\left(\mathit{r}\right(\mathit{\xi }\left)\right)=\frac{1}{2}{\langle \mathit{r}\left(\mathit{\xi }\right),\mathit{r}\left(\mathit{\xi }\right)\rangle }_{\mathbb{R}}.$$(24)

where f is the squared real norm. By taking into account the constraints (Sect. 3.2), the optimization problem finally reads:

$$\begin{array}{cc}\underset{\mathit{\xi }}{\min }& F\left(\mathit{\xi }\right)\\ \mathbf{s.t}.& {\mathit{\xi }}_{\min }\le \mathit{\xi }\le {\mathit{\xi }}_{\max }\end{array}$$(25)

where the inequality between vectors means inequality on each of their components.

3.4 Gradient computation

3.4.1 General approach

To minimize the cost function F(ξ), a gradient-based algorithm is used (e.g. quasi-Newton algorithm), which requires the computation of the gradient ∇_ξF. The Full Waveform Inversion method uses the knowledge of the forward problem to compute efficiently the gradient of the cost function. The first step is to express the explicit relation between this gradient and the derivative of the acoustic fields. As this relation involves complex operators, the complex derivative must be carefully defined. The general principles, presented in detail by Barucq et al. [32] are recalled in Appendix A. In particular, as $\mathit{\xi }\in {\mathbb{R}}^{N_{\xi }}$ and $F\left(\mathit{\xi }\right)\in \mathbb{R}$, the following relation can be shown [32] (see Appendix B):

$$\frac{\mathrm{\partial }}{\mathrm{\partial }{\xi }_j}F\left(\mathbf{\xi }\right)=2\mathfrak{R}e\left(\frac{\mathrm{d}f}{\mathrm{d}\mathbf{r}}\left(\mathbf{r}\right(\mathbf{\xi }\left)\right)\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}Z}\left(\mathbf{\xi }\right)\frac{\mathrm{d}Z}{\mathrm{d}\mathbf{U}}\left(\mathbf{\xi }\right)\frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }{\xi }_j}\left(\mathbf{\xi }\right)\right).$$(26)

The Fréchet derivative of the squared norm with respect to the residual $\frac{\mathrm{d}f}{\mathrm{d}\mathbf{r}}$ is equal to r^T, the transpose of r, and the Fréchet derivative of the impedance with respect to the acoustic field $\frac{\mathrm{d}Z}{\mathrm{d}\mathbf{U}}$ is equal to $\mathcal{R}$. The expression of the gradient is then:

$$\frac{\mathrm{\partial }}{\mathrm{\partial }{\xi }_j}F\left(\mathbf{\xi }\right)=2\mathfrak{R}e\left(\mathbf{r}\left(\mathbf{\xi }{)}^T\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}Z}\right(\mathbf{\xi }\left)\mathcal{R}\frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }{\xi }_j}\right(\mathbf{\xi })\right).$$(27)

The derivative of the residual with respect to the impedance depends on the choice of the observable. This relation is obtained from its definition equation (23) and the definition of the complex derivative (see Appendix A):

$$\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}Z}\left(\mathbf{\xi }\right)=\frac{1}{2}\left(\begin{array}{c}\left(\frac{\mathrm{d}\varphi }{\mathrm{d}Z}+\frac{\mathrm{d}\bar{\varphi }}{\mathrm{d}Z}\right)\left(Z\left(\mathbf{\xi }\right)\right)\\ -j\left(\frac{\mathrm{d}\varphi }{\mathrm{d}Z}-\frac{\mathrm{d}\bar{\varphi }}{\mathrm{d}Z}\right)\left(Z\left(\mathbf{\xi }\right)\right)\end{array}\right),$$(28)

where $\bar{\varphi }$ is the complex conjugate of ϕ. The necessity for ϕ to be ${\mathcal{C}}^1$ with respect to Z appears here clearly. Finally, the gradient of the cost function is now explicitly linked to the derivative of the acoustic fields $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}{\xi }_j}$ by equation (27).

The formal differentiation of the forward problem equation (15) with respect to any design parameter ξ_j,

$$\frac{\mathrm{d}}{\mathrm{d}{\xi }_j}\left(\mathit{A}\left(\mathbf{\xi }\right)U\left(\mathbf{\xi }\right)=\mathit{E}\right)\Rightarrow \mathit{A}\left(\mathbf{\xi }\right)\frac{\mathrm{d}U}{\mathrm{d}{\xi }_j}\left(\mathbf{\xi }\right)=-\frac{\mathrm{d}\mathit{A}}{\mathrm{d}{\xi }_j}\left(\mathbf{\xi }\right)U\left(\mathbf{\xi }\right)$$(29)

makes it explicit that this derivative can be obtained by solving a linear problem similar to equation (15) with a different source term, which could also be justified mathematically at the continuous level following a similar approach to [33]. The inverse of the matrix A being already computed (e.g., LU factorized) for the forward problem, the additional numerical cost for solving this system is negligible. The differentiation of the matrix A equation (15) follows from the differentiation of the coefficients of equations (2), (4), (8)–(10) and (12) with respect to the design parameters, along with the assembly procedure linked with the finite element method.

This approach requires solving the linear system (29) once for each design parameter. In the FWI, generally used for geophysical problems in which the number of design parameters is huge, the adjoint-state method is preferred [34], which only requires solving one additional linear system. However the benefits of this method are less important for musical instruments which involve only a few dozen design parameters. Moreover, unlike the adjoint-state method, the method presented here enables the Jacobian of the residual to be computed:

$${\nabla }_{\mathbf{\xi }}\mathbf{r}\left(\mathbf{\xi }\right)=\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}Z}\left(\mathbf{\xi }\right)\mathcal{R}\frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }{\xi }_j}\left(\mathbf{\xi }\right).$$(30)

The knowledge of this Jacobian allows the simple computation of an approximation of the Hessian: $H\approx {\nabla }_{\mathbf{\xi }}{\mathbf{r}}^T{\nabla }_{\mathbf{\xi }}\mathbf{r}$ such as used in optimization algorithms specifically suited for least-square problems (e.g. Gauss–Newton or Levenberg–Marquardt) [35], Chap. 10). In this study the optimizations are performed using the trust region reflective algorithm implemented in the Scipy python module, which takes advantage of this approximation.

3.4.2 Specific cases

This formulation is versatile with respect to the choice of the observable ϕ(Z) which only needs to be ${\mathcal{C}}^1$ with respect to Z (Eq. (28)). For example, if the reflection coefficient R_refl is chosen as the observable:

$$\varphi \left(Z\right)=R_{\mathrm{refl}}\left(Z\right):=\frac{Z-1}{Z+1}$$(31)

and defining $Z=Z_r+j{Z}_i$ with $Z_r,Z_i\in \mathbb{R}$ the real and imaginary components of Z, the derivative of ϕ is (c.f. Appendix A):

$$\frac{\mathrm{d\varphi }}{\mathrm{d}Z}=\frac{1}{2}\left(\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }{Z}_r}-j\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }{Z}_i}\right)=\frac{1}{2}\frac{1}{(Z+1{)}^2}\left(2-j\left(2j\right)\right)=\frac{2}{(Z+1{)}^2}.$$(32)

A similar computation for $\bar{\varphi }$, the conjugate of the observable, gives:

$$\{\begin{array}{ll}\frac{\mathrm{d}\varphi }{\mathrm{d}Z}=\frac{2}{(Z+1{)}^2}& \\ \frac{\mathrm{d}\bar{\varphi }}{\mathrm{d}Z}=0.& \end{array}$$(33)

This computation completes the gradient computation of the cost function with respect to the design parameters. The results for other choices of observables are given in Table 4.

4 Reconstruction strategy

The formulation of the optimization problem presented here allows the construction of the cost function to be adjusted to the targeted application. It enables a choice of the observable function ϕ but also of the considered frequencies and fingerings included in the cost function. In this study, the reconstruction of the elementary instrument with four side holes is divided into several steps corresponding to the adjustment of different design parameters. It is possible to have a different definition of the cost function at each step. The three choices (observable, frequency range and fingerings) are interlaced, each choice impacting the others. For clarity, they are here treated separately.

4.1 Choosing an observable

Following the gradient computation developed in Section 3.4, the observable must be differentiable with respect to the input impedance over the entire design domain. Ideally, an observable giving a convex function over the entire design domain (one minimum and no plateau) ensures a quick convergence of the algorithm towards a unique set of design parameters. In practice, the observable giving the widest basin of attraction and the smoothest evolution with respect to the different design variables is sought.

The observables presented in Table 4 are tested, in addition to the modulus and the unwrapped phase of the reflection coefficient. Their evolution along two dimensions of the design space are depicted in Figure 4, all other design parameters being set to their direct measured value. The chosen frequency axis ([0.1, 4] kHz with a 2 Hz step) is discussed in more detail in Section 4.2.2.

Figure 4 Costs evolution for different observables computed for frequencies in [0.1, 4] kHz with a 2 Hz step, for the fingering “oxxx” with respect to: (a) the location and, (b) the radius of the first hole. The observables are: the complex value (solid lines), the modulus (dotted lines) and the phase (dashed lines) of the impedance (gray) and the reflection coefficient (black).

Both for the location of Hole 1 (Fig. 4a) and to a lesser extent its radius (Fig. 4b), the complex impedance and its modulus show numerous fluctuations leading to many local minima in the design space. They will prevent the algorithm for converging to the correct geometry and must be excluded. This observation already made by [6], can result from the repetition of the impedance shape along the frequency axis which locally favors the concordance of two resonances even if the others do not concur. Fluctuations are also observed for the modulus of the reflection coefficient.

The unwrapped phase seems ideal for assessing the location of Hole 1 (Fig. 4a). The associated cost functions have a smooth evolution with a clear unique basin of attraction. However jumps appear for the radius of Hole 1 (Fig. 4b). This arises from the fact that the reflection coefficient crosses the origin of the complex plane for some specific geometric configurations leading to a locally non differentiable cost function. This aspect has been discussed in more detail in [17] along with a solution which requires derivation and integration with respect to the frequency axis, which is incompatible with the currently proposed formulation.

Finally the complex reflection coefficient and the impedance phase both appear acceptable in terms of convexity with the reflection coefficient being slightly preferential (Fig. 4). Similar observations can be made for all design variables taken in isolation and for different frequency ranges. The reflection coefficient defined in equation (31) is therefore chosen as the observable in this study.

4.2 Main bore geometry and hole locations

Instead of reconstructing all the design parameters together, a strategy is sought to treat the design parameters independently. As a first step, a strategy designed to recover the main bore geometry and the hole locations is investigated by using different boundary conditions given by the different fingerings.

4.2.1 Fingering choice

A design parameter can be easily adjusted in isolation, when the observable obtained with one fingering is specifically sensitive to this parameter only. To quantify this aspect, the sensitivity of each fingering with respect to the design parameters ξ_i is computed:

$${\sigma }_i\left(\mathbf{\xi }\right)=\frac{\left|\right|{\mathrm{\partial }}_{{\xi }_i}\varphi \left(Z\right(\mathbf{\omega },\mathbf{\xi }\left)\right)|{|}_{L2,\mathbf{\omega }}}{\left|\right|\varphi \left(Z\right(\mathbf{\omega },\mathbf{\xi }\left)\right)|{|}_{L2,\mathbf{\omega }}}.$$(34)

The sensitivity with respect to the locations of the holes and to the main bore length for each fingering is plotted in Figure 5. The computation is carried out for the direct measured geometry. This sensitivity map remains globally similar regardless of the value of the design parameters, as long as the topology of the instrument is unchanged (order of the holes, etc.). All the fingerings are sensitive to the main bore length and, to a lesser extent, to the location of the open holes, particularly the first ones. The locations of the upstream closed holes have a small effect on the observable (1 order of magnitude smaller). The downstream closed holes have negligible effect.

Figure 5 Sensitivity of the reflection coefficient w.r.t the hole locations/bore length design parameters (Eq. (34)) for frequencies in [0.1, 4] kHz with a 2 Hz step. Open (resp. closed) holes are indicated with circles (resp. crosses).

Following these observations, the locations of the holes have a good degree of independence, if only fingerings with one open hole are used (xxxo, xxox, xoxx and oxxx). The location of each hole can be recovered by performing independent optimizations on the corresponding fingering (hole 4 from xxxo, hole 3 from xxox, etc.), after having first recovered the main bore length from the all closed fingering.

This process can only provide a rough estimation of the geometry because the locations of the holes cannot be estimated independently of their radius and chimney height, which in turn cannot be recovered using a unique fingering (see Sect. 4.3). With respect to the main bore, its length is strongly inter-linked with its output radius (conicity). Hopefully, this radius and length can be estimated together with sufficient accuracy from the all closed fingering (“xxxx”).

4.2.2 Frequency range

As the observable and the fingering are chosen to recover each hole location and the main bore geometry, it is interesting to investigate how the choice of the frequency range can facilitate the convergence of the optimization algorithm.

Figure 6 compares the evolution of the cost function computed with $\varphi \left(Z\right)=R_{\mathrm{refl}}$ on different sets of frequencies with respect to the main bore length (Fig. 6a) and the location of the first hole (Fig. 6b), the other parameters being fixed at their directly measured values.

Figure 6 Cost evolution along one dimension of the design space for three frequency ranges: [0.1, 0.5] kHz and [0.1, 4] kHz with a 100 Hz step, and [0.1, 4] kHz with a 2 Hz step. (a) Total length Lmb, fingering: xxxx’’. (b): Location of hole 1 Xh1, fingering “oxxx’’.

The first observation is that the frequency step has little influence on the cost value. It can be increased from 2 Hz to 100 Hz with negligible effect on the cost function (Fig. 6). Contrary to the usual assumption, the physical information is not concentrated at the resonance, rather it influences the acoustic response at any frequency. The number and the distribution of the observed frequencies only need to be chosen to avoid confusion between geometries which could have similar responses at one given frequency and to reduce the effect of the measurement noise on the reconstructed geometry. The choice of a 100 Hz step appears here a good compromise between computation cost and the accuracy of the reconstruction.

Secondly, in both cases it appears that reducing the frequency range to [100 Hz, 500 Hz] increases the basin of attraction. This can arise from the fact that, at high frequencies, two resonances can accidentally match for a given geometry leading to a local minimum; this will not be the case at low frequency, especially for frequencies lower than the first resonance. It is noticeable that ultimately, only 5 frequency values are necessary to obtain a satisfactory smooth cost-function. The conclusion of the two previous sections on the observable choice (Sect. 4.1) and on the fingering choice (Sect. 4.2.1) are still valid with this frequency range.

4.3 Hole dimensions

The conclusions of Section 4.2 can not be applied to the radii and the chimney height s of the holes. There is no sense in treating these two dimensions separately, as they have a similar effect on the impedance at low frequencies. Figure 7 shows the influence of the frequency range on the cost function evolution with respect to the first hole radius (a similar effect can be observed with the chimney height). Here again, refining the frequency step does not modify the cost function evolution, a 100 Hz step is sufficient. But this time, the low frequencies are not sufficient: the basin of attraction is so large that the minimum almost disappears. The inclusion of higher frequencies allows the minimum to be refined.

Figure 7 Cost evolution with respect to the radius of Hole 1 on the fingering “oxxx”, with three frequency ranges: [0.1, 0.5] kHz and [0.1, 4] kHz with a 100 Hz step, and [0.1, 4] kHz with a 2 Hz step.

Simultaneously recovering both the radius of the hole and its chimney height is a challenge. To clearly distinguish their effects, it is necessary to reach frequencies such that $k{H}_{\mathrm{hi}}>1$, which corresponding to f > 27 kHz for 2 mm for the presented instrument. In this range, the modeling is not valid anymore and we are not able to measure the impedance. However, these two parameters differ also by the fact that only the radius appears in the expression for the hole radiation. Accounting for several boundary conditions through the different fingerings will help in recovering their values. This is illustrated by the observation of the cost function in the design space plane (R _h1, H _h1) shown in Figure 8. When only one fingering (oxxx) is considered, a valley without minimum appears, whereas when fingerings including different open-closed states for Hole 1 are considered, a clear minimum appears. To conclude, recovering the dimensions of the holes necessitates consideration over both a wide frequency range and several fingerings.

Figure 8 Cost corresponding to the reflection coefficient, with respect to the first hole radius Rh1 and chimney height Hh1, for frequencies in [0.1, 4] kHz with a 0.1 kHz step. (a) The fingering for which only the first hole is open (oxxx). (b) Five fingerings (xxxx, xxxo, xxox, xoxx, oxxx)

4.4 Results

Following the analysis of Sections 4.2 and 4.3, a two-step strategy is proposed to reconstruct the instrument geometry.

1. A rough estimation is first performed by recovering only the main bore geometry and the hole locations over low frequencies only ([100, 500] Hz, step 100 Hz).

2. Then, a finer reconstruction is performed by including all design parameters over a wider frequency range ([100, 4000] Hz, step 100 Hz).

The rough estimation is divided into five optimization problems. The main bore geometry (L_mb and R_mb(L_mb)) is optimized from the all closed fingering (xxxx). Then each hole location is recovered, from the last one to the first one, by observing the fingering for which only the hole of interest is open (“xxxo” for Hole 4, etc.). Each optimization converges in about 8 iterations and the total rough estimation is computed in about 5 s on a personal computer. After this step, the hole locations and total length deviate by about 5 mm from their measured values and the radius R _mb deviates by about 0.05 mm. The refined reconstruction is then performed by considering the five fingerings (“xxxx”, “xxxo”, “xxox”, “xoxx”, “oxxx”), and activating the adjustment of the 14 design parameters. The optimization algorithm converges after about 15 iterations computed in about 1 min on a personal computer. To evaluate the robustness of the process, the total reconstruction is performed on three sets of data obtained from three repeated measurements of the impedance of each fingering.

The geometry obtained for the first set of data is compared to the initial state and the direct measurements in Figure 9. The resulting deviations for the three sets of data are presented in more detail in Figure 10. The direct measurement uncertainties (see Tab. 1) are indicated as dotted lines. The reconstruction uncertainties are estimated by using the sensitivity of the observable σ_i defined in equation (34) for the resulting geometry and the impedance measurement accuracy estimated in Section 2.6 ($\Delta \varphi \left(Z\right)=5\%$):

$$\Delta {\xi }_i=\frac{\Delta \varphi \left(Z\right)}{{\sigma }_i}.$$(35)

Figure 9 Comparison between the geometry obtained by direct measurement (wide black line), the initial state (thin black line) and the reconstructed geometry from the first set of data (grey). For readability the vertical and horizontal axis have different scales.

Figure 10 Deviation between manually measured and reconstructed geometries from the three sets of data. The crosses represent the results for the mixed radiation condition (holes flanged, main bore unflanged) and the triangle the results with only unflanged radiation for the first set of data. The direct measurement uncertainties (see Tab. 1) are indicated as dotted lines. (a) Hole locations and main bore length. (b) Hole radii and main bore radius. (c) Chimney lengths.

Since the sensitivity depends on the fingering, the uncertainty Δξ_i is computed for each of them and the smallest value is depicted.

The geometric values obtained by reconstruction are in good agreement with the those obtained by direct measurements. The deviations on the hole locations and the main bore length are within 0.5 mm (Fig. 10a), within 0.1 mm for the radii (Fig. 10b) and within 0.3 mm for the chimney lengths (Fig. 10c). In terms of the uncertainties and the accuracy needed to rebuild an instrument, these results are satisfactory.

The variation in the results between the sets of data is more pronounced for the hole dimensions (Figs. 10b and 10c). Hole 1 is poorly reconstructed with the third set of data. This could be due to a measurement issue such as a leak. The holes are closed with a piece of Teflon covered with putty. Even when applied carefully, the putty can slightly penetrate the chimney, resulting in a different effective length between open and closed states and distorting the reconstruction. As the tube is made of brass, its temperature is very sensitive to handling, which can also alter the measurements.

In general, the hole radii appear to be slightly underestimated (Fig. 10b). This bias could come from the model assumptions such as the radiation condition (Sect. 2.4). To observe this influence, a reconstruction has been performed on the first set of data with unflanged radiation for all the openings (black triangles on Fig. 10b). This primarily reduces the radii and increases the chimney heights. It therefore deteriorates the results. It is possible to use more refined radiation models [36]. As the main pipe external radius is 4 mm, cylindrical flange or finite flange conditions should result here in an intermediary condition in between unflanged and infinite flanged conditions at low frequency. However, these radiation conditions strongly depend on frequency and could eventually improve the reconstruction.

5 Conclusion

In this study, the Full Waveform Inversion technique is used for the bore reconstruction of woodwind-like musical instruments. This methodology exploits the fact that the gradient of the cost function can be computed by solving the same linear system which arises from the finite element spatial discretization of the direct problem of wave propagation. This gradient computation is more accurate and more efficient than finite-difference. It therefore improves the efficiency of gradient-based algorithms performing the inversion. The approach developed here allows the computation of the gradient of the cost function with respect to any design variable and whatever the geometry of the instrument (variable cross section of the main bore, side hole). The only limitation is that the observable used in the cost function must be differentiable with respect to the acoustic fields (pressure and/or flow). This is one of the main developments of this work with respect to previous studies.

The analysis of the cost functions constructed from different observables reveals that the reflection coefficient is a better candidate than the input impedance for reconstructing the bore profile. Only a few frequencies are necessary to correctly estimate the geometric values. To recover the hole locations, the exclusion of high frequencies increases the basin of attraction of the cost function. To recover the hole dimensions (radius and chimney height) it is necessary to include the high frequencies and to account for several boundary conditions (different fingerings). Following these considerations, the use of local optimization algorithms to perform the reconstruction appears suitable. The phase of the reflection coefficient appears to be a good candidate to further improve the algorithm’s convergence, in spite of the fact that a loss of regularity occurs due to its wrapping. This aspect should be considered further in future work.

This method is applied successfully to reconstruct an elementary woodwind-like instrument from measured data. For this purpose, a strategy including step by step inversion on specific fingerings is devised. This is a promising step from the perspective of using this method in the context of manufacture. It could allow estimation of the internal geometry of historical instruments or aid the design of new instruments. The remaining limitations preventing the use of this method on real instruments is the complexity of the main bore shape, limited here to a simple conical part. One interesting approach could be to combine pulse reflectometry on the all closed fingering to have a first estimation of the main bore shape, before performing the FWI strategy proposed in this study. Reconstruction of a real instrument also involves the reconstruction of much more than four side holes. This should not complicate the rough estimation of their location as each hole is treated independently. However all holes should be treated together to estimate the chimneys dimensions. Increasing both the number of design variables and the number of fingerings included in the cost function will significantly increase the computation time. It could also impede the optimization process. The presence of keywork over the holes will also complicate the reconstruction. If it can not be avoided (by removing the key for the impedance measurement), their effect should be added to the radiation expressions for the side holes and the related design variables should be added to the optimization problem. Since the effect of keywork is close to a change in the chimney geometry [36], the resulting inverse problem should be harder to solve.

In seismology, it is well identified that the FWI gives better results if more measurement points are used instead of increasing the number of frequencies. Since our impedance sensor uses five microphones, another interesting approach could be to reconstruct the instrument directly from the five pressure signals instead of the experimentally derived input impedance.

6 Data availability

All tools used in this study have been implemented in the open-source software OpenWInD [20]. The scripts and the data are available under GPLv3 licence at https://gitlab.inria.fr/openwind/openwind.

Appendix A

Complex derivative of complex function

If $z\in \mathbb{C}$, the derivative operator $\frac{\mathrm{d}}{\mathrm{d}z}$ needs to be precisely defined. This definition follows the work of Barucq et al. [32] inspired by the work of Brandwood [37]. Each application $g:\mathbb{C}\to \mathbb{C}$ is interpreted as the application $g:\mathbb{C}\times \mathbb{C}\to \mathbb{C}$ of the complex number z and its conjugate $\bar{z}$, with $z=z_r+j{z}_i$ and $\bar{z}=z_r-j{z}_i$. The variables z and $\bar{z}$ are then considered as independent variables. The application $f:\mathbb{R}\times \mathbb{R}\to \mathbb{C}$ is defined from g such as $f(z_r,z_i)=g(z,\bar{z})$. The complex derivatives are defined from the real derivative $\mathrm{\partial }/\mathrm{\partial }{z}_r$ and $\mathrm{\partial }/\mathrm{\partial }{z}_i$:

$$\frac{\mathrm{\partial }g}{\mathrm{\partial }z}=\frac{1}{2}\left(\frac{\mathrm{\partial }f}{\mathrm{\partial }{z}_r}-j\frac{\mathrm{\partial }f}{\mathrm{\partial }{z}_i}\right),$$(36)

$$\frac{\mathrm{\partial }g}{\mathrm{\partial }\bar{z}}=\frac{1}{2}\left(\frac{\mathrm{\partial }f}{\mathrm{\partial }x}+j\frac{\mathrm{\partial }f}{\mathrm{\partial }{z}_i}\right).$$(37)

Remark: In the general case where $g(z,\bar{z})\in \mathbb{C}$, it can be noticed that:

$$\overline{\frac{\mathrm{\partial }g}{\mathrm{\partial }\stackrel{̅}{z}}}=\overline{\frac{1}{2}\left(\frac{\mathrm{\partial }f}{\mathrm{\partial }{z}_r}+j\frac{\mathrm{\partial }f}{\mathrm{\partial }{z}_i}\right)}=\frac{1}{2}\left(\frac{\mathrm{\partial }\bar{f}}{\mathrm{\partial }{z}_r}-j\frac{\mathrm{\partial }\bar{f}}{\mathrm{\partial }{z}_i}\right)=\frac{\mathrm{\partial }\bar{g}}{\mathrm{\partial }z}$$(38)

which is an extension of the corollary A.1 given in the appendix of the paper [32].

Appendix B

Complex derivative of real function

Let

$$\begin{array}{llll}h:& R& \to & R\\ & x& \to & h\left(x\right)\equiv g\left(z\right(x),\bar{z}(x\left)\right)\end{array},$$(39)

where $g:{\mathbb{C}}^2\to \mathbb{R}$ and $z:\mathbb{R}\to \mathbb{C}$. By applying the definition of equation (37) we can show that (see [32]):

$$\frac{\mathrm{d}h}{\mathrm{d}x}=\mathfrak{R}e\left(\frac{\mathrm{\partial }g}{\mathrm{\partial }z}\frac{\mathrm{\partial }z}{\mathrm{\partial }x}+\frac{\mathrm{\partial }g}{\mathrm{\partial }\bar{z}}\frac{\mathrm{\partial }\bar{z}}{\mathrm{\partial }x}\right)=\mathfrak{R}e\left(\frac{\mathrm{\partial }g}{\mathrm{\partial }z}\frac{\mathrm{\partial }z}{\mathrm{\partial }x}\right)+\mathfrak{R}e\left(\frac{\mathrm{\partial }g}{\mathrm{\partial }\bar{z}}\frac{\mathrm{\partial }\bar{z}}{\mathrm{\partial }x}\right).$$(40)

Applying the result of equation (38) twice, we can remark that:

$$\mathfrak{R}e\left(\frac{\mathrm{\partial }g}{\mathrm{\partial }\bar{z}}\frac{\mathrm{\partial }\bar{z}}{\mathrm{\partial }x}\right)=\mathfrak{R}e\left(\overline{\frac{\mathrm{\partial }\stackrel{̅}{g}}{\mathrm{\partial }z}}\frac{\mathrm{\partial }\bar{z}}{\mathrm{\partial }x}\right)=\mathfrak{R}e\left(\frac{\mathrm{\partial }\bar{g}}{\mathrm{\partial }z}\overline{\frac{\mathrm{\partial }\stackrel{̅}{z}}{\mathrm{\partial }x}}\right)$$(41)

$$=\mathfrak{R}e\left(\frac{\mathrm{\partial }\bar{g}}{\mathrm{\partial }z}\frac{\mathrm{\partial }z}{\mathrm{\partial }\bar{x}}\right).$$(42)

In this specific case, where $x\in \mathbb{R}$ and $g(z,\bar{z})\in \mathbb{R}$ we finally get:

$$\frac{\mathrm{d}h}{\mathrm{d}x}=2\mathfrak{R}e\left(\frac{\mathrm{\partial }g}{\mathrm{\partial }z}\frac{\mathrm{\partial }z}{\mathrm{\partial }x}\right).$$(43)

Acknowledgments

This research has been supported by the Nouvelle-Aquitaine regional council (Cofimus project).

References

All Tables

All Figures

	Figure 1 Sketch of the cylindrical tube with four side holes to be recovered. It is described by 15 geometric parameters (three for the main bore profile and three for each of the four holes).
In the text

	Figure 2 The network representation of the studied resonator is composed of 9 pipes and 10 connection conditions.
In the text

	Figure 3 Sketch of a T-junction with the adopted convention for the flow orientation.
In the text

	Figure 4 Costs evolution for different observables computed for frequencies in [0.1, 4] kHz with a 2 Hz step, for the fingering “oxxx” with respect to: (a) the location and, (b) the radius of the first hole. The observables are: the complex value (solid lines), the modulus (dotted lines) and the phase (dashed lines) of the impedance (gray) and the reflection coefficient (black).
In the text

	Figure 5 Sensitivity of the reflection coefficient w.r.t the hole locations/bore length design parameters (Eq. (34)) for frequencies in [0.1, 4] kHz with a 2 Hz step. Open (resp. closed) holes are indicated with circles (resp. crosses).
In the text

	Figure 6 Cost evolution along one dimension of the design space for three frequency ranges: [0.1, 0.5] kHz and [0.1, 4] kHz with a 100 Hz step, and [0.1, 4] kHz with a 2 Hz step. (a) Total length Lmb, fingering: xxxx’’. (b): Location of hole 1 Xh1, fingering “oxxx’’.
In the text

	Figure 7 Cost evolution with respect to the radius of Hole 1 on the fingering “oxxx”, with three frequency ranges: [0.1, 0.5] kHz and [0.1, 4] kHz with a 100 Hz step, and [0.1, 4] kHz with a 2 Hz step.
In the text

	Figure 8 Cost corresponding to the reflection coefficient, with respect to the first hole radius Rh1 and chimney height Hh1, for frequencies in [0.1, 4] kHz with a 0.1 kHz step. (a) The fingering for which only the first hole is open (oxxx). (b) Five fingerings (xxxx, xxxo, xxox, xoxx, oxxx)
In the text

	Figure 9 Comparison between the geometry obtained by direct measurement (wide black line), the initial state (thin black line) and the reconstructed geometry from the first set of data (grey). For readability the vertical and horizontal axis have different scales.
In the text

Figure 10 Deviation between manually measured and reconstructed geometries from the three sets of data. The crosses represent the results for the mixed radiation condition (holes flanged, main bore unflanged) and the triangle the results with only unflanged radiation for the first set of data. The direct measurement uncertainties (see Tab. 1) are indicated as dotted lines. (a) Hole locations and main bore length. (b) Hole radii and main bore radius. (c) Chimney lengths.

In the text