Issue 
Acta Acust.
Volume 5, 2021



Article Number  53  
Number of page(s)  13  
Section  Musical Acoustics  
DOI  https://doi.org/10.1051/aacus/2021046  
Published online  03 December 2021 
Technical & Applied Article
Peakpicking identification technique for modal expansion of input impedance of brass instruments
Laboratoire d’Acoustique de l’Université du Mans, CNRS – UMR 6613, avenue Olivier Messiaen, 72085 Le Mans Cedex 9, France
^{*} Corresponding author: frederic.ablitzer@univlemans.fr
Received:
12
July
2021
Accepted:
2
November
2021
The paper presents a method to obtain the modal expansion of the measured input impedance of a brass instrument. The method operates as a peakpicking procedure, which makes it particularly intuitive for users who are not experts in modal analysis. To bypass the limitation of usual peakpicking approaches, which are valid only for well separated resonances, the present method is based on a semilocal optimization problem. It consists in adjusting the frequency and damping of one mode at a time while taking into account the presence of all other modes in the basis. The practical application of the method involves four elementary actions, which can be chained in different ways to progressively approximate a measured input impedance. This procedure is illustrated through the approximation of the input impedance of a bass trombone. The supervised nature of the method allows the user to favour modes that have a physical meaning, i.e. that can be associated with a resonance peak. A single spurious mode can however be deliberately introduced to approximate the input impedance curve beyond the last visible peak. The method applies directly to the frequencydomain data provided by an impedance sensor and does not require any preprocessing. Nevertheless, it is fairly robust to noisy data. Since the method allows a reconstruction of the input impedance using either complex modes or real modes, results obtained with each approximation are critically compared.
Key words: Modal analysis / Peakpicking / Input impedance / Brass instrument
© F. Ablitzer, Published by EDP Sciences, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The input impedance of a wind instrument can be mathematically represented using a modal expansion [1]. By carrying out modal identification from a measured input impedance, it is therefore possible to model the resonator of an existing instrument by a set of ordinary differential equations which can be directly plugged into numerical methods for analysing the emergence of selfsustained oscillations, such as timedomain simulations [2], linear stability analysis [3–5] or numerical continuation [6–8]. This approach is facilitated by the fact that wellestablished experimental devices and procedures exist for measuring the input impedance of wind instruments [9, 10]. In contrast, there is no consensus on how to obtain the modal parameters from the measured input impedance. Various homemade algorithms are proposed to identify modal parameters, some of which are inspired from methods typically used for the experimental modal analysis of vibrating structures (for an overview of such methods, the reader may refer for example to [11]).
Modal identification techniques are commonly classified into two categories: single degree of freedom (SDOF) and multiple degrees of freedom (MDOF) methods. The main assumption behind SDOF methods is that the input impedance (or more generally any frequency response function) can be well approximated near a resonance by the theoretical response of a single degree of freedom system. Wellknown SDOF methods are the −3 dB method, based on estimating the halfpower bandwidth of a resonance peak, or the circle fitting method, which exploits the fact that the frequency response of a SDOF system around its resonance depicts a circle in the Nyquist plane [12]. Such methods are also commonly described as peakpicking methods, since they use the data around a resonance peak to obtain the parameters of the corresponding mode. They are simple to implement and easy to operate by users who are not experts in modal analysis. However, they perform poorly if the underlying assumption is not met, i.e. when neighbouring modes significantly contribute to the response around a resonance peak (this phenomenon is described as modal overlap). The procedure described in [13] to obtain modal parameters from a measured input impedance belongs to this category of modal identification techniques. The author emphasizes the interest of this approach, which requires no curve fitting, while noting that it is vulnerable to modal overlap in high frequency, leading to a significant deviation of the approximation from the original input impedance.
In contrast, MDOF methods attempt to identify several modal components simultaneously within a given frequency range. They are more robust to the effect of modal overlap, but generally less readily available and less simple to implement. MDOF identification techniques are based on a mathematical representation of the system in terms of poles (carrying information on natural frequencies and modal damping ratios) and residues (which represents the amplitude of modal components). A common approach in the practical application of MDOF techniques is a twostep process [11]: in a first step, the poles of the system are estimated; in a second step, these poles are used to obtain the residues. This approach has been adopted in most of the abovementionned studies (e.g. [4, 6, 7]). In a recent paper dedicated to modal analysis of wind instruments [14], Taillard et al. also propose a twostep method, which combines the Least Squares Complex Exponential (LSCE) technique to identify the poles with a fitting procedure to obtain the residues. The proposed method requires a number of preprocessing steps to be applicable to experimental data, possibly noisy and available over a limited frequency range. The fitting procedure leads to an excellent approximation of the measured input impedance, at the price of including more poles in the basis than the number of resonance peaks visible on the curve.
In most contexts where experimental modal analysis is carried out, the determination of the order of the model (i.e., the number of poles to consider) is a key issue. Most commercially available modal analysis software packages leave this decision to the user, by providing a widely used representation named “stabilization diagram” [11]. This decisionsupport tool is very useful but can be difficult to exploit for users who are not experts in modal analysis. A radically opposite approach is to estimate the optimal order without involving the user. This approach can be described as unsupervised modal analysis, which is still considered as a challenging problem [15], especially when the modal overlap is high in the frequency range of interest. As an attempt to overcome this issue, high resolution modal analysis techniques have been developed. One approach sometimes used in musical acoustics is the ESPRIT technique (Estimation of Signal Parameters via Rotational Invariance Techniques) [16], which enables the identification of poles in situations where the modal overlap is strong. This technique is often combined with the ESTER (Estimation Error) criterion [17] to determine the optimal order of the modal model. Although mainly applied to vibrating structures (e.g. piano soundboard [18], guitar [19]), this approach has recently been applied to brass instruments [6], leading to a satisfactory reconstruction of the measured input impedance with a limited number of poles (each one being associated with a resonance peak). In a methodology developed by Maestre et al. for modelling a measured input impedance as a recursive digital filter [20], another strategy is proposed to set the order of the model: the lowfrequency region is described with as many poles as resonance peaks, while the highfrequency region is described with a set of poles uniformly distributed on a logarithmic frequency axis (the positions of all poles are then optimized iteratively).
Alternatively to the classification into SDOF/MDOF methods, modal identification techniques can be classified into frequency domain and time domain methods, depending on the form in which the measured data must be presented to enable the identification of poles. For instance, the ESPRIT technique is a time domain method, in the sense that it exploits the impulse response of the system, expected to be a sum of exponentially decaying sinusoids. Such a method can be applied quite straightforwardly, or with little preprocessing, when the response of the system is obtained with impulse excitation (typically, impact hammer excitation) [16]. On the other hand, when the system is excited using broadband excitation (shaker excitation for vibrating structures, impedance sensor for acoustic resonators), its response is obtained in the form of a frequency response function (drivingpoint mobility, input impedance). As a result, a preliminary step in the application of a time domain modal identification technique is the conversion of the input impedance into an impulse response, which typically involves several preprocessing operations (extrapolation of data outside the measurement frequency range, filtering, etc.) [6].
The present paper proposes an alternative modal identification technique, which aims at combining the simplicity of SDOF methods and their intuitive interpretation (“each resonance peak corresponds to a mode”) with the robustness of MDOF methods, where several modes are assumed to contribute to the response near each resonance. The proposed method operates as a peakpicking procedure, in the sense that the identification of each mode consists in “capturing” a resonance peak visible on the measured input impedance.
The main features of the proposed approach are listed below:
it is supervised, in the sense that it requires actions from the user;
it is iterative: the procedure involves a stepbystep enrichment of the modal basis and progressive refinement of the reconstructed input impedance;
it is a frequencydomain method, which directly exploits the measured input impedance;
it requires no preprocessing of the measured data (such as filtering, low frequency extrapolation, etc.) and is fairly robust to noisy data;
it allows an approximation of the input impedance using either complex or real modes.
Regarding the last point, the literature review shows that the representation using real modes seems less popular than that using complex modes to approximate a measured input impedance [4, 6, 14]. The motivation to prefer complex modes is however rarely stated explicitly, except when the dynamics of the resonator is studied in a theoretical framework, where the complex nature of modes results from modeling assumptions (in particular concerning losses) [21]. To the best of our knowledge, a comparison between the two approaches has never been shown in a case where it is desired to reconstruct the input impedance up to high frequencies (i.e. including the last visible, although “weak”, resonance peaks). Particular attention is paid to this issue in the paper.
The principle of the identification technique is presented in Section 2. Section 3 gives indications on how to implement the method. A typical way to operate the method is illustrated through the reconstruction of the input impedance of a brass instrument in Section 4. Key points in the manner of formulating the modal identification problem are also analysed and discussed in this section.
2 Principle of the method
2.1 Modal expansion of the input impedance
The normalized input impedance Z(ω) of an acoustic resonator can be represented as a sum of complex modes [1],
$$Z\left(\omega \right)=\sum _{n=1}^{N}\mathrm{}\left(\frac{{c}_{n}}{\mathrm{j\omega}{s}_{n}}+\frac{{c}_{n}^{\mathrm{*}}}{\mathrm{j\omega}{s}_{n}^{\mathrm{*}}}\right).$$(1)
Considering this form of modal expansion of the input impedance, the identification problem consists in finding two complex unknowns for each mode, namely the pole s_{n} and the associated residue c_{n}. The poles s_{n} are related to natural frequencies ω_{n} and modal damping ratios ξ_{n} as [22]:
$${s}_{n}={\xi}_{n}{\omega}_{n}+j{\omega}_{n}\sqrt{1{\xi}_{n}^{2}}.$$(2)
The modal decomposition in equation (1) can be rewritten as:
$$Z\left(\omega \right)=\sum _{n=1}^{N}\mathrm{}\frac{\mathrm{j\omega}{A}_{n}+{B}_{n}}{{\omega}_{n}^{2}{\omega}^{2}+j2{\xi}_{n}{\omega}_{n}\omega},$$(3)where
$$\{\begin{array}{l}{A}_{n}=2\mathcal{R}\mathrm{e}\left({c}_{n}\right)\\ {B}_{n}=2\mathcal{R}\mathrm{e}\left({c}_{n}{s}_{n}^{\mathrm{*}}\right),\end{array}$$(4)are realvalued amplitudes, from which the complex residues c_{n} can reciprocally be calculated if necessary as:
$$\{\begin{array}{l}\mathcal{R}e\left({c}_{n}\right)=\frac{{A}_{n}}{2}\\ \mathrm{Im}\left({c}_{n}\right)=\frac{1}{2}\frac{{\xi}_{n}{\omega}_{n}{A}_{n}{B}_{n}}{{\omega}_{n}\sqrt{1{\xi}_{n}^{2}}}.\\ \end{array}$$(5)
The identification problem therefore consists in finding four real unknowns for each mode (ω_{n}, ξ_{n}, A_{n}, B_{n}) such that equation (3) provides a fair approximation of the measured input impedance.
Alternatively, a representation of the input impedance as a sum of real modes can be considered. The corresponding modal decomposition is [1]
$$Z\left(\omega \right)=\sum _{n=1}^{N}\mathrm{}\frac{\mathrm{j\omega}{A}_{n}}{{\omega}_{n}^{2}{\omega}^{2}+j2{\xi}_{n}{\omega}_{n}\omega}.$$(6)
In this case, the identification problem consists in finding three real unknowns for each mode (ω_{n}, ξ_{n}, A_{n}).
2.2 Estimation of modal amplitudes
Assuming first that the number N of poles and their values (ω_{n}, ξ_{n}) are known, the problem of finding the real constants A_{n} and B_{n} (if complex modes are considered) reduces to solving a system of linear equations, which can be written in matrix form as:
$$\mathbf{M}\left(\mathbf{\nu},{\omega}_{1},{\xi}_{1},\dots ,{\omega}_{N},{\xi}_{N}\right)\mathbf{x}=\mathbf{d}\left(\mathbf{\nu}\right),$$(7)where
$$\mathbf{M}(\mathbf{\nu},{\omega}_{1},{\xi}_{1},\dots ,{\omega}_{N},{\xi}_{N})=\left[\begin{array}{ll}\mathcal{R}e\left({\mathbf{M}}_{A}\right)& \mathcal{R}e\left({\mathbf{M}}_{B}\right)\\ \mathcal{I}m\left({\mathbf{M}}_{A}\right)& \mathcal{I}m\left({\mathbf{M}}_{B}\right)\end{array}\right]$$(8)with
$$\{\begin{array}{l}{\mathbf{M}}_{A}[i,n]=\frac{j\mathbf{\nu}\left[i\right]}{{\omega}_{n}^{2}\mathbf{\nu}[i{]}^{2}+j2{\xi}_{n}{\omega}_{n}\mathbf{\nu}[i]}\\ {\mathbf{M}}_{B}[i,n]=\frac{1}{{\omega}_{n}^{2}\mathbf{\nu}[i{]}^{2}+j2{\xi}_{n}{\omega}_{n}\mathbf{\nu}[i]},\end{array}$$(9)
$$\mathbf{d}\left(\mathbf{\nu}\right)=\left[\begin{array}{l}\mathcal{R}e\left(Z\left(\mathbf{\nu}\right)\right)\\ \mathcal{I}m\left(Z\left(\mathbf{\nu}\right)\right)\end{array}\right]$$(10)and
$$\mathbf{x}=\left[\begin{array}{l}{A}_{1}\\ \vdots \\ {A}_{N}\\ {B}_{1}\\ \vdots \\ {B}_{N}\end{array}\right].$$(11)
Explicitly partitioning M and d into real and imaginary parts of complex quantities guaranties that the amplitudes A_{n} and B_{n} obtained from a least squares solution are real valued and therefore respect the proper form of modal decomposition (Eq. (3)). In equations (8) and (10), ν represents a column vector containing the angular frequencies considered in the identification of amplitude coefficients. The choice of an appropriate frequency range will be discussed in Section 3.4. Matrix M is written as a block matrix for compactness in equation (8). As indicated by equation (9), the rows of each block correspond to the frequencies contained in ν, whereas the columns correspond to the modal components. If K frequency lines are considered to identify the amplitude of N modes, matrix M has dimension (2K, 2N) (with K ≫ N in general). Equation (7) is therefore solved in a least squares sense as
$$\widehat{\mathbf{x}}(\mathbf{\nu},{\omega}_{1},{\xi}_{1},\dots ,{\omega}_{N},{\xi}_{N})=\underset{\mathbf{x}}{\mathrm{argmin}}\left\right\mathbf{M}(\mathbf{\nu},{\omega}_{1},{\xi}_{1},\dots ,{\omega}_{N},{\xi}_{N})\mathbf{x}\mathbf{d}\left(\mathbf{\nu}\right){}^{2}$$(12)to obtain the modal amplitudes (the hat (̂) over a quantity represents its estimate).
If it is chosen to approximate the input impedance using real modes, the same equation has to be solved, yet with different expressions of matrix M and vector x,
$$\mathbf{M}(\mathbf{\nu},{\omega}_{1},{\xi}_{1},\dots ,{\omega}_{N},{\xi}_{N})=\left[\begin{array}{l}\mathcal{R}e\left({\mathbf{M}}_{A}\right)\\ \mathcal{I}m\left({\mathbf{M}}_{A}\right)\end{array}\right]$$(13)and
$$\mathbf{x}=\left[\begin{array}{l}{A}_{1}\\ \vdots \\ {A}_{N}\end{array}\right],$$(14)where M_{A} has the same expression as previously given in equation (9).
It should be noted that the values of A_{n} obtained with matrix M defined by equation (13) can differ from those obtained with matrix M defined by equation (8), even if the same poles are considered. In other words, considering real modes in the estimation of amplitudes with equation (12) is not equivalent to setting each B_{n} to 0 in a solution obtained by considering complex modes.
2.3 Estimation of poles
In contrast to obtaining the amplitudes if the poles are known, the determination of the poles themselves cannot be reduced to solving a linear system. In absolute terms, the number N of poles necessary to obtain a satisfactory approximation of the input impedance should be considered as an unknown of the problem. However, we first address the problem in which it is assumed to be known, or set by the user. Under this assumption, the identification problem can be written as:
$$({\widehat{\omega}}_{1},{\widehat{\xi}}_{1},\dots ,{\widehat{\omega}}_{N},{\widehat{\xi}}_{N})=\underset{({\omega}_{1},{\xi}_{1},\dots ,{\omega}_{N},{\xi}_{N})}{\mathrm{argmin}}J(\mathbf{\nu},{\omega}_{1},{\xi}_{1},\dots ,{\omega}_{N},{\xi}_{N}),$$(15)where the cost function
$$J(\mathbf{\nu},{\omega}_{1},{\xi}_{1},\dots ,{\omega}_{N},{\xi}_{N})=\left\right\mathbf{M}(\mathbf{\nu},{\omega}_{1},{\xi}_{1},\dots ,{\omega}_{N},{\xi}_{N})\widehat{\mathbf{x}}(\mathbf{\nu},{\omega}_{1},{\xi}_{1},\dots ,{\omega}_{N},{\xi}_{N})\mathbf{d}\left(\mathbf{\nu}\right){}^{2}$$(16)corresponds to the residual sum of squares associated with the estimation of amplitudes (Eq. (12)). Equation (15) represents a global optimization problem, which consists in simultaneously finding the values of all poles and the corresponding amplitudes (which are estimated at each evaluation of the cost function). Although it is expected to provide the best estimate of modal parameters, it has several drawbacks:
the computational effort to solve the problem is high, since the minimum of a function of 2N variables has to be found,
there is no guarantee that the optimization process will not end up stuck in a local minimum,
the cost function cannot be displayed with usual modes of representations (such as a curve, a surface or volume slices for functions of 1, 2 or 3 variables, respectively), which reduces the possibilities of diagnosis if difficulties are encountered during the identification process.
To overcome these issues, we introduce an auxiliary problem, expressed mathematically as:
$$({\widehat{\omega}}_{n},{\widehat{\xi}}_{n})=\underset{({\omega}_{n},{\xi}_{n})}{\mathrm{argmin}}J(\mathbf{\nu},{\widehat{\omega}}_{1},{\widehat{\xi}}_{1},\dots ,{\omega}_{n},{\xi}_{n},\dots ,{\widehat{\omega}}_{N},{\widehat{\xi}}_{N}).$$(17)
This problem can be referred to as a semilocal optimization problem. It consists in finding the value of one pole, specified by rank n in equation (17), while the value of all other poles is fixed. In contrast to the previous problem, it involves much less computational effort and allows a visualization of the cost function, since the latter is a function of only two variables. Although this optimization problem is local with respect to a pole, it should be emphasized that the calculation of amplitudes at each evaluation of the cost function is done for all modes considered in the sum. In other words, during the optimization process of adjusting one pole through equation (17), the value of all other poles are left unchanged, but the associated amplitudes can vary. This way of formulating the problem is motivated by the substantial modal overlap that often characterizes the last resonance peaks in the input impedance of a brass instrument. In this frequency region, the amplitude to be assigned to a modal component may depend strongly on the characteristics of the neighbouring modes (amplitudes and poles). Adjusting the amplitude of all modal components during the optimization of a pole is therefore a desirable feature.
The general principle of the proposed identification procedure is to replace the global optimization problem (Eq. (15)) by a succession of semilocal optimization problems (Eq. (17)), where the modal components needed to reconstruct the input impedance are introduced one after the other, in a usercontrolled manner. As will be shown later, this approach involves an iterative refinement of the reconstructed input impedance. At this stage, it should be pointed out that there is no argument for asserting that a succession of semilocal problems converge towards the global solution of Eq. (15). Particular attention will be paid to this issue in Section 4.2.
3 Practical implementation of the method
The proposed peakpicking identification technique is based on four elementary actions:
adding a mode (A),
removing a mode (R),
refining the identification (F),
changing the fitting bandwidth (B).
This section describes the purpose of each action and how to implement it. A typical way to chain them to obtain the modal model of a measured input impedance will be illustrated in Section 4.2.
3.1 Adding a mode
Adding a mode is the main action of the method. It consists in capturing a resonance peak visible on the input impedance by introducing a new pole in the basis and solving equation (17) in order to adjust its natural frequency ω_{n} and damping ratio ξ_{n}. The initial guess for the natural frequency can conveniently be provided by clicking a resonance peak visible on the amplitude curve, as most scientific computing environments provide a user interaction function to do so (e.g. ginput in Matlab or matplotlib.pyplot.ginput with Python). The initial guess for the damping ratio can be considered as a hidden parameter, from the user point of view. A constant value of 5% regardless of mode rank is a reasonable choice for a brass instrument [23].
In practice, finding the minimum of the cost function involved in equation (17) can be done using the NelderMead simplex algorithm, readily available in common scientific computing environments (e.g. function fminsearch in Matlab or scipy.optimize.minimize with Python). Solving equation (7) in a least squares sense within the cost function is computationally inexpensive and straightforward, since it is a linear system of equations. It can be done for instance using mldivide in Matlab or linalg.lstsq with Python.
Thanks to the semilocal formulation of the optimization problem, the parameters of the current mode (frequency, damping, amplitude) are found taking into account all other modes already present in the basis, if any. Conversely, the amplitudes associated with all other modes are simultaneously adapted to the presence of the new one.
3.2 Removing a mode
Unsurprisingly, this action consists in removing one of the poles from the modal basis. Immediately after this operation, the amplitudes associated with the remaining poles are then updated using the least squares solution of equation (7).
This action is useful when the attempt to capture a resonance peak fails and produces a spurious pole [24], i.e. with a natural frequency very far from the initial guess, often associated with an unusually high damping ratio. When this occurs, however, it may be useful to keep the spurious pole during one or more iterations before removing it, since its presence in the basis can help to capture the desired peaks, as will be illustrated later.
3.3 Refining the identification
This action consists in solving equation (17) successively for each mode already present in the basis. It is motivated by the fact that the optimum associated with one pole is likely to change after introduction of additional poles into the basis, due to the semilocal nature of the optimization problem. The refinement action is then useful to update the value of each pole taking into account the current value of other poles in the basis. In a very concrete way, refining the identification is done by solving equation (17) within a loop from the first to the last modal index (the modes being sorted in ascending order of natural frequency, regardless of the order in which the resonance peaks were captured). Updating the poles in ascending order of natural frequency during the refinement process has been found to provide more robustness (see the companion website^{1} for an illustration).
It must be stressed that the optimum found when updating each pole through equation (17) is a temporary optimum: as soon as the value of the next pole in the list will be modified, itself through an optimization process, the cost function associated with the present pole will have changed. Performing the refinement action several times and observing the convergence is a means of assessing the mutual sensitivity between the estimated poles. A way to do so will be proposed in Section 4.2.
3.4 Changing the fitting bandwidth
As evoked in Section 2.2, the estimation of modal amplitudes through equation (7) involves only a subset of frequency lines among all those included in the measurement data. The considered frequency lines are represented by vector ν. The aim is to minimize the difference between the measured input impedance and its modal approximation only near the resonance peaks that the user decides to capture.
Practically, each time a mode is added to capture a resonance peak, a frequency range centred around the peak is added to the vector ν. More precisely, the interval associated with each mode is [f_{n,0} – Δf/2, f_{n,0} + Δf/2], where f_{n,0} is the initial guess for the natural frequency, the value of which is set when adding a mode. The total frequency range considered in the estimation of amplitudes is the union of all intervals (there are as many as there are modes in the basis). It should be noted that the total frequency range can be discontinuous. This happens if Δf is less than the spacing between two adjacent resonance peaks.
In the following, the parameter Δf controlling the width of local intervals is named “fitting bandwidth”. This parameter is generally set to a low value at the beginning of the identification procedure; once a sufficient number of modal components have been obtained, the value can be increased such that a continuous frequency range is considered during the last iterations. What is meant by a low value depends on the resonator, and primarily on the length of its bore, which controls the spacing between the resonance peaks. In the practical implementation of the method, it is advisable to highlight graphically the current frequency range considered in the identification, as a visual aid for setting an appropriate value of Δf.
4 Application
4.1 Approximation of the input impedance of a bass trombone
The instrument considered to illustrate the application of the method is a Courtois bass trombone [25]. Its input impedance measured with the slide in first position is shown in Figure 1. The magnitude of the input impedance exhibits a series of “strong” resonances (1–9), followed by “weak” resonances (10–12). Each of these resonances, except the first one, allows the musician to play a note with a pitch close to the resonance frequency (such a note is referred to as a “natural note”) [26]. The 12th resonance peak corresponds to a F5 note. Although rarely used in the repertoire, this very high note is playable by skilled players and can therefore be considered within the pitch range of the instrument. Beyond these peaks, two “very weak” resonances can be distinguished (13 and 14).
Figure 1 Measured input impedance (black thick line) of a bass trombone and its approximation using 15 complex modes (blue line). Individual modal contributions are shown as thin gray lines. The difference between the measured input impedance Z_{mes} and its approximation Z (red line) is expressed as Z  Z_{mes} in (a) and arg(Z)  arg(Z_{mes}) in (b). 
Figure 1 shows an approximation of the input impedance using 15 complex modes, obtained with the proposed method. The corresponding modal parameters (natural frequencies, damping ratios and amplitudes) are listed in Table 1. To each of the resonance peaks just described corresponds a mode, clearly observable in the magnitude plot (see Fig. 1a). A 15th mode, also visible through its individual contribution, has been chosen to be included into the basis. Although this mode cannot definitely be associated with a resonance peak, it contributes to achieve a better agreement between the modal model and the experimental data in the upper frequency range. As can be seen, the proposed method leads to a satisfactory reconstruction of the input impedance. Importantly, it is obtained with a limited number of poles in the basis, with the intention of being able to associate almost each mode with a resonance peak.
Regarding the phase plot (see Fig. 1b), it can be noticed that the reconstructed input impedance satisfies the passivity condition (phase comprised in the [−π/2, π/2] interval, such that $\mathcal{R}e\left(Z\right)\ge 0$), whereas some individual modal components do not. In particular, the phase associated with mode 11 to mode 15, which correspond to weak resonances, systematically exceed π/2 over a wide frequency range. Such behaviour is associated with negative values of A_{n} and/or B_{n} in the least squares estimation of amplitudes (see Tab. 1). Similarly, the phase of the reconstructed input impedance converges towards zero at 0 Hz (in the present solution, $Z\left(0\right)={\sum}_{n=1}^{N}\mathrm{}{B}_{n}/{\omega}_{n}^{2}$ is a positive real value), whereas the phase of some modal components do not (those associated with negative values of B_{n}).
4.2 Iterative reconstruction of the input impedance
The reconstructed input impedance shown in Figure 1 is the result of a sequence of 23 elementary actions. A way to represent this sequence is proposed in Figure 2a, which indicates the successive actions together with the agreement between the measured input impedance and its approximation. The chosen indicator is the value of the cost function (Eq. (16)) normalized by the number K of frequency lines considered,
$$\mathrm{res}=\frac{J(\mathbf{\nu},{\widehat{\omega}}_{1},{\widehat{\xi}}_{1},\dots ,{\widehat{\omega}}_{N},{\widehat{\xi}}_{N})}{K}.$$(18)
Figure 2 (a) Evolution of the residual associated with the least squares solution of equation (7) throughout the iterative reconstruction of the input impedance shown in Figure 1. Each iteration corresponds to the execution of one action by the user, among the four possible ones. Each action is labeled by a letter (see the text for the description of each action). (b) Evolution of poles, represented by their natural frequencies f_{n} and damping ratios ξ_{n}, throughout the iterative approximation of the input impedance. The color indicates the progression of this process, using the same colormap as in (a). The locations of poles at iteration 23, corresponding to the result shown in Figure 1, are highlighted (blue circles). Vertical dotted lines indicate a harmonic series based on the second resonance frequency. 
To facilitate the interpretation of this graph, the actions are explicitly described hereafter:
iterations 1–12 (A): the first 12 resonance peaks are successfully captured by adding modes;
iter. 13 (A): a mode is added to try capturing the 13th resonance peak, which results in a spurious pole (f_{n} ≈ 180 Hz, ξ_{n} ≈ 67%);
iter. 14 (R): the spurious pole is removed;
iter. 15 (F): refinement is performed;
iter. 16–17 (A): the 13th and 14th resonance peaks are now successfully captured by adding modes;
iter. 18 (B): the fitting bandwidth Δf is enlarged from 40 Hz to 80 Hz;
iter. 19 (F): refinement is performed;
iter. 20 (A): a mode is added to try capturing a hypothetical 15th resonance peak, which results in a spurious pole (f_{n} ≈ 1380 Hz, ξ_{n} ≈ 10%);
iter. 21 (A): keeping the spurious pole, a 15th resonance is now successfully captured by adding a mode;
iter. 22 (R): the spurious pole is removed;
iter. 23 (F): refinement is performed; in practice, the process can be stopped after this iteration, since a satisfactory approximation of the input impedance has been obtained (see Fig. 1);
iter. 24–70 (F): numerous additional refinement iterations are performed, for illustrative purposes.
In general, strong resonances are captured on the first try, with a very low risk of obtaining a spurious pole. In the case illustrated here, the first spurious pole occurs when trying to capture the 13th peak, which corresponds to a very weak resonance. Such failures can be avoided (or repaired) either by performing refinement more frequently, or by voluntarily keeping a spurious pole in the basis until the targeted resonances are successfully captured. These two scenarios have been illustrated here (the corresponding sequences are iterations 13–17 and 20–22, respectively).
In the practical application of the method, the graphical representation proposed in Figure 2a can be used as a decision aid to find out when it is no longer useful to refine the approximation, after the desired number of resonances has been captured. In the present case, the residual only slightly diminishes after iteration 23, which means that the approximation is not significantly improved by a succession of additional refinement steps. It must be emphasized that the sequence for achieving a satisfactory reconstruction is not unique: another sequence of actions may lead to a similar end result. A comparison between 4 different sequences is proposed on the companion website (see Footnote ^{1}).
The evolution of poles during the iterative reconstruction of the input impedance is shown in Figure 2b. The locations of the first nine poles remain almost unchanged throughout the process. This is consistent with the fact that the corresponding resonances are strong, so that the individual contribution of one mode is dominant close to each peak (see Fig. 1a). In contrast, the locations of poles associated with weaker resonances (10–14) significantly vary during the process, which indicates that their optimum value with respect to the semilocal optimization problem (Eq. (17)) is affected by other poles in the basis. This is consistent with what is observed in Figure 1a for resonances 10–14: the amplitude of the reconstructed input impedance in the vicinity of each peak differs from the amplitude of the corresponding modal component, as well as the frequency at which the maximum is reached. This observation reveals why usual peakpicking methods, based on a SDOF assumption, would not be suitable to obtain a modal model that includes all visible peaks. The interdependency between neighbouring poles is especially visible on the last three poles (13–15). While these poles follows the same trend at iteration 23 (proximity to the harmonic series and similar values of damping ratio), performing numerous further refinement iterations makes the 15th pole slowly move away to a higher frequency and much higher damping ratio, which in turn modifies the locations of the 13th and 14th poles.
Interestingly, the final locations of poles 13 and 14 seem to be physically meaningful. Their natural frequencies indeed fit well within the quasi harmonic series formed by the resonances of lower rank. Likewise, the values of damping ratio are similar to those obtained for previous resonances. This observation indicates that intentionally including a last spurious pole in the modal basis can help to identify modes with the highest possible physical meaning (this would be useful, for example, if they were to be compared with the modes obtained from a numerical model of the resonator). From this point of view, the spurious pole plays a role similar to that of the upper residual correction often considered in experimental modal analysis of structures to represent the effect of outofband modes [11]. On the contrary, in the absence of a spurious pole, the optimization algorithm adjusts the last poles such as to compensate the mismatch related to the truncation of the modal series.
4.3 Mutual sensitivity between poles
To further illustrate the principle of the semilocal optimization problem, the cost function involved in equation (17) is examined for two cases.
First, the cost function associated with the 9th mode, corresponding to the last strong resonance, is considered. Figure 3a shows the cost function immediately after this mode has been introduced (iter. 9), whereas Figure 3b shows the cost function after all other modes have been added (iter. 23). The cost function remains almost the same between these two distant iterations, which indicates that the optimum value of this pole is not sensitive to that of other poles (it should not be concluded, however, that the associated amplitude remains the same between iterations 9 and 23). This behaviour is typical for cost functions associated with strong resonances.
Figure 3 Cost function associated with the 9th mode at (a) iteration 9 (b) iteration 23. The scale ranges from black (lowest value) to white (highest value). The white circle indicates the initial guess, the green asterisks indicates the minimum found by the optimization algorithm. 
Then, the cost function associated with the 13th mode, corresponding to a very weak resonance, is analysed. Figure 4a shows the cost function immediately after this mode has been successfully introduced (iter. 16). The location of the minimum is not obvious. Instead, the large dark region indicates that many values of the pole tend to provide a good agreement between the modal model and the experimental data. More specifically, the shape of the cost function can be interpreted as follows: the optimization algorithm might “hesitate” between capturing the resonance peak (region slightly below the initial guess in Fig. 4a), or producing a spurious mode with a higher frequency and excessive damping ratio (region on the right of the initial guess in Fig. 4a). Figure 4b shows the same cost function after the next mode, corresponding to the 14th resonance peak, has been added (iter. 17). The presence of a new mode changes the shape of the cost function. The location of the minimum appears more clearly. It is also better insulated from regions likely to attract the optimization algorithm towards values far from the initial guess. This example shows that the presence of at least one higher mode acts as a favourable constraint to keep a pole within an admissible range of values during the optimization process, such that it can be associated with a resonance peak rather than transforming into a spurious pole.
Figure 4 Cost function associated with the 13th mode at (a) iteration 16 (b) iteration 17. The scale ranges from black (lowest value) to white (highest value). The white circle indicates the initial guess, the green asterisks indicates the minimum found by the optimization algorithm. 
4.4 Choosing between complex modes or real modes
Until now, the proposed identification method has been illustrated using a representation of the input impedance as a sum of complex modes. Alternatively, a sum of real modes can be considered without changing the way the method is applied in practice. It is even possible to switch from one representation to the the other at any step of the procedure, simply with a recalculation of the least squares solution of equation (7) using one matrix M or the other (Eq. (8) for complex modes, Eq. (13) for real modes). In practice, it is useful to perform at least one refinement action after this, since the values of poles allowing an optimal reconstruction of the input impedance may slightly differ between the two representations.
The input impedance of the same instrument is now approximated using a sum of real modes. Figure 5 compares the results obtained with each approximation. The modal parameters associated with the real modes approximation are listed in Table 2. A first look at magnitude and phase curves reveals that both complex modes and real modes provide a fair approximation of the measured input impedance. Regarding the residual, it cannot be stated that one representation is clearly surpassing the other. Interestingly, the approximation with real modes yields a slightly lower residual in the upper frequency range, in the vicinity of weak resonance peaks, which is rather counterintuitive. Since complex modes offer more freedom in the curve fitting problem expressed by equation (7) (two unknown coefficients for the amplitude of each complex mode, against one for each real mode), one could expect that they provide a better approximation than real modes in a frequency region where the modal overlap is high. The present result indicates on the contrary that the improvement occurs mainly in the low frequency region. In light of these observations, the added value of using complex modes is actually low for the instrument considered here.
Figure 5 Comparison between the reconstructed input impedance of the bass trombone using 15 complex modes (blue continuous line) or 15 real modes (orange dashed line). The difference between the measured input impedance and each approximation is shown in (a) using the same line style. The individual modal contributions shown as thin gray lines correspond to the real modes approximation. 
Regarding the phase plot (see Fig. 5b), it can be observed that the phase of each component is comprised in the [−π/2, π/2] interval. With the real modes approximation, the passivity condition is therefore verified not only for the reconstructed input impedance, but also for each individual modal contribution (including the spurious mode), as opposed to the result obtained with complex modes (see Fig. 1b). Another difference between the two approximations, which goes unnoticed on the magnitude plot with the chosen scale but is observable on the phase plot, concerns the behaviour of the reconstructed input impedance when the frequency approaches zero. Whereas the representation with real modes always leads to a zero value (as seen in Eq. (3) by setting B_{n} = 0 and ω = 0), the representation with complex modes permits a nonzero value at ω = 0 ($Z\left(0\right)={\sum}_{n=1}^{N}\mathrm{}{B}_{n}/{\omega}_{n}^{2}$). Although it is possible to give it a physical meaning, it is hazardous to trust the value at ω = 0 resulting from any curvefitting procedure (including the present one), since measurement data at very low frequencies is typically either missing or highly corrupted with noise. For this reason, it seems advisable to systematically adopt a real mode representation to approximate the input impedance of brass instruments. Another argument in favour of this recommendation is given in the next section.
4.5 Robustness to noisy measurements
In order to assess the sensitivity of the method to measurement uncertainties, the reconstruction of the input impedance of another instrument is carried out. The instrument is a contrabass tuba (16 ft instrument), which has a conical bore longer than the bass trombone (9 ft instrument) previously considered. As a result, its input impedance, shown in Figure 6, exhibits a series of resonance peaks which is shifted towards lower frequencies compared to the previous instrument, by a factor about 1.78. In particular, the first two peaks are located in a frequency range where the signaltonoise ratio is especially poor (see Figs. 6c and 6d), due to the lower level the excitation source of the impedance sensor can provide.
Figure 6 Reconstructed input impedance of a contrabass tuba using 15 complex modes (blue continuous line) and 15 real modes (orange dashed line). The arrow in (d) indicates a frequency region where the reconstructed input impedance using complex modes violates the passivity condition. 
The results of modal identifications using 15 complex modes or 15 real modes are also shown in Figure 6. The first 14 modes correspond to individual resonance peaks. A 15th spurious mode has been introduced to mimic the effect of several highly overlapped modes beyond 500 Hz, which are difficult to capture as individual components without artificially increasing the number of poles in the basis. In spite of the poor signaltonoise ratio at low frequencies, the peakpicking procedure performs successfully without having to preprocess the measurement data, regardless of the chosen approximation (complex or real modes).
In the case of complex modes, however, although the reconstructed input impedance is overall very satisfactory, closer inspection at low frequencies reveals phase values slightly outside the [−π/2, π/2] interval (see Fig. 6d), i.e. such that $\mathcal{R}e\left(Z\right)<0$, which is an undesired result. For instance, this modal model could not be safely used to study selfsustained oscillation regimes, as it violates the passivity condition. The inappropriate behavior of the identification with complex modes can be explained by the greater freedom of this modal model, which makes it more likely to fit the noise in the data. Although passivity of the reconstructed input impedance could be enforced through preprocessing of the measured data, which is the approach proposed in [14], it is advisable to take advantage of the regularizing effect of a real modes approximation when the measurement is very noisy.
5 Conclusion
An identification technique for the modal analysis of input impedance of brass instruments was developed. Rather than seeking automatic extraction of modal parameters from the measured data, the proposed method operates as a peakpicking procedure: the reconstruction of the input impedance is done progressively by successively capturing visible resonance peaks on the experimental curve, in a supervised manner.
The method was shown to achieve a satisfying reconstruction of the input impedance. In particular, it is able to capture the last, weak resonance peaks typically exhibited by brass instruments. Being able to associate each mode with a resonance peak can be particularly interesting for the subsequent analysis of the instrument. Once the modal basis has been obtained, it is possible for example to change the amplitude associated with a given resonance peak to study its influence of the selfsustained oscillation regimes.
The robustness of the proposed identification technique to measurement noise was demonstrated. The choice between a representation with complex modes or real modes was also discussed. In light of the results obtained, it seems that an identification based on real modes is to be preferred for brass instruments. Although the representation using complex modes has a more general scope, it does not bring a significantly better agreement between the measured input impedance and its modal approximation. Moreover, it is likely to bring undesirable effects. In particular, a case where the input impedance reconstructed with complex modes does not respect the passivity condition has been highlighted.
Up to now, the method presented in the paper has been applied to various brass instruments. Some cases are shown for illustration purposes on a companion website (Footnote ^{1}). The extent to which the method is also applicable to instruments with open toneholes, for which the input impedance exhibits less regularly spaced resonance peaks, remains to be examined. While remaining cautious, it is hoped that the proposed method will also prove useful for the modal analysis of woodwind instruments.
Conflict of interest
Author declared no conflict of interests.
Acknowledgments
The author would like to thank Joël Gilbert, JeanPierre Dalmont and Sylvain Maugeais for helpful comments on a draft of this manuscript.
References
 A. Chaigne, J. Kergomard: Acoustics of musical instruments. Springer, 2016. [CrossRef] [Google Scholar]
 F. Silva, C. Vergez, P. Guillemain, J. Kergomard, V. Debut: Moreesc: a framework for the simulation and analysis of sound production in reed and brass instruments. Acta Acustica united with Acustica 100, 1 (2014) 126–138. [CrossRef] [Google Scholar]
 J.S. Cullen, J. Gilbert, D.M. Campbell: Brass instruments: linear stability analysis and experiments with an artificial mouth. Acta Acustica united with Acustica 86, 4 (2000) 704–724. [Google Scholar]
 L. Velut, C. Vergez, J. Gilbert, M. Djahanbani: How well can linear stability analysis predict the behaviour of an outwardstriking valve brass instrument model?. Acta Acustica united with Acustica 103, 1 (2017) 132–148. [CrossRef] [Google Scholar]
 V. Debut, J. Antunes, O. Inácio: Linear modal stability analysis of bowedstrings. The Journal of the Acoustical Society of America 141, 3 (2017) 2107–2120. [CrossRef] [PubMed] [Google Scholar]
 V. Fréour, L. Guillot, H. Masuda, S. Usa, E. Tominaga, Y. Tohgi, C. Vergez, B. Cochelin: Numerical continuation of a physical model of brass instruments: Application to trumpet comparisons. The Journal of the Acoustical Society of America 148, 2 (2020) 748–758. [CrossRef] [PubMed] [Google Scholar]
 J. Gilbert, S. Maugeais, C. Vergez: Minimal blowing pressure allowing periodic oscillations in a simplified reed musical instrument model: BouasseBenade prescription assessed through numerical continuation. Acta Acustica 46 (2020) 27. [CrossRef] [EDP Sciences] [Google Scholar]
 T. Colinot, L. Guillot, C. Vergez, P. Guillemain, J.B. Doc, B. Cochelin: Influence of the “ghost reed” simplification on the bifurcation diagram of a saxophone model. Acta Acustica united with Acustica 105, 6 (2019) 1291–1294. [CrossRef] [Google Scholar]
 C.A. Macaluso, J.P. Dalmont: Trumpet with nearperfect harmonicity: Design and acoustic results. The Journal of the Acoustical Society of America 129, 1 (2011) 404–414. [CrossRef] [PubMed] [Google Scholar]
 W. Kausel: Bore reconstruction of tubular ducts from its acoustic input impedance curve. IEEE Transactions on Instrumentation and Measurement 53, 4 (2004) 1097–1105. [CrossRef] [Google Scholar]
 P. Avitabile: Modal testing: a practitioner’s guide. John Wiley & Sons, 2017. [CrossRef] [Google Scholar]
 D.J. Ewins: Modal testing: theory, practice and application. John Wiley & Sons, 2009. [Google Scholar]
 T. Colinot: Numerical simulation of woodwind dynamics: investigating nonlinear sound production behavior in saxophonelike instruments, PhD thesis. AixMarseille Université, 2020. [Google Scholar]
 P.A. Taillard, F. Silva, P. Guillemain, J. Kergomard: Modal analysis of the input impedance of wind instruments. application to the sound synthesis of a clarinet. Applied Acoustics 141 (2018) 271–280. [CrossRef] [Google Scholar]
 B. Peeters, J. Lau, J. Lanslot, H. Van der Auweraer: Automatic modal analysis – myth or reality? Sound and Vibration 42, 3 (2008) 17. [Google Scholar]
 K. Ege, X. Boutillon, B. David: Highresolution modal analysis. Journal of Sound and Vibration 325, 4–5 (2009) 852–869. [CrossRef] [Google Scholar]
 R. Badeau, B. David, G. Richard: A new perturbation analysis for signal enumeration in rotational invariance techniques. IEEE Transactions on Signal Processing 54, 2 (2006) 450–458. [CrossRef] [Google Scholar]
 K. Ege, X. Boutillon, M. Rébillat: Vibroacoustics of the piano soundboard: (non) linearity and modal properties in the lowand midfrequency ranges. Journal of Sound and Vibration 332, 5 (2013) 1288–1305. [CrossRef] [Google Scholar]
 B. Elie, F. Gautier, B. David: Macro parameters describing the mechanical behavior of classical guitars. The Journal of the Acoustical Society of America 132, 6 (2012) 4013–4024. [CrossRef] [PubMed] [Google Scholar]
 E. Maestre, G.P. Scavone, J.O. Smith, Joint modeling of impedance and radiation as a recursive parallel filter structure for efficient synthesis of wind instrument sound, in: Proceedings of the 21st International Conference on Digital Audio Effects, 2018, pp. 4–8. [Google Scholar]
 J. Kergomard, V. Debut, D. Matignon: Resonance modes in a onedimensional medium with two purely resistive boundaries: Calculation methods, orthogonality, and completeness. The Journal of the Acoustical Society of America 119, 3 (2006) 1356–1367. [CrossRef] [Google Scholar]
 M. Géradin, D.J. Rixen: Mechanical vibrations: theory and application to structural dynamics. John Wiley & Sons, 2014. [Google Scholar]
 A. MamouMani, D.B. Sharp: Evaluating the suitability of acoustical measurement techniques and psychophysical testing for studying the consistency of musical wind instrument manufacturing. Applied Acoustics 71, 7 (2010) 668–674. [CrossRef] [Google Scholar]
 I. Markovsky, J. Boets, B. Vanluyten, K. De Cock, B. De Moor: When is a pole spurious? in: International Conference on Noise and Vibration Engineering, 2006, pp. 1615–1626. [Google Scholar]
 J. Gilbert, L. Leblanc, C. Vergez: L’analyse de stabilité linéaire pour évaluer la facilité d’émission des cuivres. etude comparative de trombones ténor et basse, in: 14e Congrès Français d’Acoustique, Le Havre, France, 2018. [Google Scholar]
 M. Campbell, J. Gilbert, A. Myers: The science of brass instruments, Springer, 2021. [CrossRef] [Google Scholar]
Cite this article as: Ablitzer F. 2021. Peakpicking identification technique for modal expansion of input impedance of brass instruments. Acta Acustica, 5, 53.
All Tables
All Figures
Figure 1 Measured input impedance (black thick line) of a bass trombone and its approximation using 15 complex modes (blue line). Individual modal contributions are shown as thin gray lines. The difference between the measured input impedance Z_{mes} and its approximation Z (red line) is expressed as Z  Z_{mes} in (a) and arg(Z)  arg(Z_{mes}) in (b). 

In the text 
Figure 2 (a) Evolution of the residual associated with the least squares solution of equation (7) throughout the iterative reconstruction of the input impedance shown in Figure 1. Each iteration corresponds to the execution of one action by the user, among the four possible ones. Each action is labeled by a letter (see the text for the description of each action). (b) Evolution of poles, represented by their natural frequencies f_{n} and damping ratios ξ_{n}, throughout the iterative approximation of the input impedance. The color indicates the progression of this process, using the same colormap as in (a). The locations of poles at iteration 23, corresponding to the result shown in Figure 1, are highlighted (blue circles). Vertical dotted lines indicate a harmonic series based on the second resonance frequency. 

In the text 
Figure 3 Cost function associated with the 9th mode at (a) iteration 9 (b) iteration 23. The scale ranges from black (lowest value) to white (highest value). The white circle indicates the initial guess, the green asterisks indicates the minimum found by the optimization algorithm. 

In the text 
Figure 4 Cost function associated with the 13th mode at (a) iteration 16 (b) iteration 17. The scale ranges from black (lowest value) to white (highest value). The white circle indicates the initial guess, the green asterisks indicates the minimum found by the optimization algorithm. 

In the text 
Figure 5 Comparison between the reconstructed input impedance of the bass trombone using 15 complex modes (blue continuous line) or 15 real modes (orange dashed line). The difference between the measured input impedance and each approximation is shown in (a) using the same line style. The individual modal contributions shown as thin gray lines correspond to the real modes approximation. 

In the text 
Figure 6 Reconstructed input impedance of a contrabass tuba using 15 complex modes (blue continuous line) and 15 real modes (orange dashed line). The arrow in (d) indicates a frequency region where the reconstructed input impedance using complex modes violates the passivity condition. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.