Pierre-André Taillard

Franck Laloë

Michel Gross

Jean-Pierre Dalmont

Jean Kergomard

Statistical estimation of mechanical parameters of clarinet reeds using experimental and numerical approaches

A set of 55 clarinet reeds is observed by holography, collecting 2 series of measurements made under 2 different moisture contents, from which the resonance frequencies of the 15 first modes are deduced. A statistical analysis of the results reveals good correlations, but also significant differences between both series. Within a given series, flexural modes are not strongly correlated. A Principal Component Analysis (PCA) shows that the measurements of each series can be described with 3 factors capturing more than 90% of the variance: the first is linked with transverse modes, the second with flexural modes of high order and the third with the first flexural mode. A forth factor is necessary to take into account the individual sensitivity to moisture content. Numerical 3D simulations are conducted by Finite Element Method, based on a given reed shape and an orthotropic model. A sensitivity analysis revels that, besides the density, the theoretical frequencies depend mainly on 2 parameters: EL and GLT . An approximate analytical formula is proposed to calculate the resonance frequencies as a function of these 2 parameters. The discrepancy between the observed frequencies and those calculated with the analytical formula suggests that the elastic moduli of the measured reeds are frequency dependent. A viscoelastic model is then developed, whose parameters are computed as a linear combination from 4 orthogonal components, using a standard least squares fitting procedure and leading to an objective characterization of the material properties of the cane Arundo donax.

Introduction

Clarinettists experience every day the crucial importance of clarinet reeds for the quality of sound. Their characterization is a real challenge for musicians who wish to obtain reeds that are suited to their personal needs. The present paper address this complex field of research. Its scope is restricted to the development of an objective method for a mechanical characterization of single reeds of clarinet type. >From the shape and the resonance frequencies of each individual reed (measured with heterodyne holography), we intend to deduce the mechanical properties of the material composing it. A subsequent study should then examine how these mechanical properties are correlated with the musical properties of the reeds.

Generally, the physicist chooses a model in order to validate it by observations. In the present study, the complexity of the problematic forced us to adopt the reverse attitude: We observe the mechanical behavior of clarinet reeds with a statistically representative sample and exploit afterward the statistical results for establishing a satisfactory mechanical model designed with a minimal number of parameters.

Natural materials, as wood or cane, are often orthotropic and exhibit a different stiffness along the grain (longitudinally) as in the others directions. The problem is then obviously multidimensional. Nevertheless, reed makers classify their reeds by a single parameter: the nominal reed "strength" (also called "hardness"), in general from 1 to 5, which basically reflects the stiffness of the material (cane, Arundo donax L.), since all reeds of the same model have theoretically the same shape. The method of measurement is generally not publicized by manufacturers, but this "strength" is probably related to the static Young modulus in the longitudinal direction E L .

"Static" (i.e. low frequency) measurements of the elastic parameters of cane are available in the literature, for instance Spatz et al. [START_REF] Spatz | Biomechanics of the giant reed Arundo donax[END_REF]. A viscoelastic behavior has been reported in experimental situations (see e.g. Marandas et al. [START_REF] Marandas | Caractérisation mécanique des anches simples d'instruments à vent[END_REF], Ollivier [3] or Dalmont et al. [START_REF] Dalmont | Nonlinear characteristics of single-reed instruments: Quasistatic volume flow and reed opening measurements[END_REF]) and this fact seems generally well accepted in wood sciences and biomechanics (for instance Speck et al. [START_REF] Speck | Mechanical Properties of the Rhizome of Arundo donax L[END_REF][START_REF] Speck | Damped oscillations of the giant reed Arundo donax (Poaceae)[END_REF]). Marandas et al. proposed a viscoplastic model of the wet reed. Viscoelastic behavior for cane was already demonstrated by Chevaux [START_REF] Ph | Améliorations de la durée de vie des anches doubles en roseau pour instruments à vents[END_REF], Obataya et al. [START_REF] Obataya | Physical properties of cane used for clarinet reed[END_REF][START_REF] Obataya | The effects of water soluble extractives on the acoustic properties of reed (Arundo donax L.)[END_REF][START_REF] Obataya | Mechanical relaxation processes due to sugars in cane (Arundo donax L.)[END_REF][START_REF] Obataya | Acoustic properties of a reed (Arundo donax L.) used for the vibrating plate of a clarinet[END_REF] and Lord [START_REF] Lord | Viscoelasticity of the giant reed material Arundo donax[END_REF]. These authors study only the viscoelasticity of the longitudinal Young modulus E L , leaving aside the case of the shear modulus in the longitudinal/tangential plane G LT . Furthermore, they give no really representative statistics about the variability of the measured parameters.

The observation of mechanical resonance frequencies can be achieved by different methods. The methods used by Chevaux, Obataya and Lord are destructive for the reed, which cannot be used for further musical tests. On the contrary, holography is a convenient non-destructive method, the reed being excited by a loudspeaker. For instance Pinard et al. [START_REF] Pinard | Musical quality assessment of clarinet reeds using optical holography[END_REF] measured with this method the frequency of the 4 lowest resonances and focused their attention on the musical properties of the reeds.

The digital Fresnel holography method was used by Picart et al. [START_REF] Picart | Tracking high amplitude auto-oscillations with digital Fresnel holograms[END_REF][START_REF] Picart | Study of the Mechanical Behaviour of a Clarinet Reed Under Forced and Auto-oscillations With Digital Fresnel Holography[END_REF] and Mounier et al. [START_REF] Mounier | Investigation of clarinet reed auto-oscillations with digital Fresnel holography[END_REF] to measure high amplitude motion of a reed blown by an artificial mouth. Guimezanes [START_REF] Guimezanes | Etude Expérimentale et Numérique de l'Anche de Clarinette[END_REF] used a scanning vibrometer.

Recent technological developments provide very efficient and convenient measurements with holography, without having to manually identify the modes of resonance and to be satisfied with a single picture of their vibration: in a few minutes hundreds of holograms are acquired showing the response of a reed for many frequencies. The temperature and the moisture content can be considered as constant during a measurement series 1 . The Sideband Digital Holography technique provides additional facilities (see 2.1.1).

Different authors (among them Casadonte [START_REF] Casadonte | The perfect clarinet reed? Vibrational modes of realistic clarinet reeds[END_REF][START_REF] Casadonte | The clarinet reed: an introduction to its biology, chemistry and physics[END_REF], Facchinetti et al. [START_REF] Facchinetti | Application of modal analysis and synthesis of reed and pipe to numerical simulations of a clarinet[END_REF][START_REF] Facchinetti | Numerical and experimental modal analysis of the reed and pipe of a clarinet[END_REF] and Guimezanes [START_REF] Guimezanes | Etude Expérimentale et Numérique de l'Anche de Clarinette[END_REF]) modeled the clarinet reed by Finite Elements Method (FEM) and computed the first few eigenmodes. They chose appropriated values of the elastic parameters in the literature, ignoring however viscoelastic behavior. The goodness of fit between observations and model was of secondary importance, except for Guimezanes. This latter author built a 2-D elastic model of the reed with longitudinally varying parameters. He fitted his model quite adequately with his observations (only 5 resonances were measured), but the fitted parameters seem not really plausible physically. His model didn't respect the assumption of a radial monotonically decrease of stiffness from the outer side to the inner side of the cane. Under such conditions, the frequency of the first resonance would increase in comparison to homogeneous material, and not decreased, as observed experimentally.

In Section 2 the measurement method is presented. The experimental setup is described in Section 2.1 and the method for observing resonance frequencies is detailed in 2.2. The results for 55 reeds are given in Section 3 (statistics, Principal Component Analysis (PCA) [START_REF] Jolliffe | Principal component analysis[END_REF]).

In Section 4, the development and the selection of a satisfactory mechanical model with minimal structure is described. First, a numerical analysis of the resonance frequencies of a reed assumed to be perfectly elastic is done by Finite Element Method (FEM), and a metamodel computing the resonance frequencies from elastic parameters is given in Section 4.3. This allows solving the inverse problem in a fast way. However, because the elastic model is not very satisfactory, viscoelasticity has to be introduced and some parameters are added to the model in Section 4.4. The viscoelastic model has however too many degrees of freedom, according to PCA. Consequently, the viscoelastic parameters of the model are assumed to be correlated and PCA indicates that these parameters can be probably reconstructed from 4 orthogonal components, as a linear combination, by multiple regression (Section 5). The relationships between the components and the viscoelastic parameters is given, and finally the resulting values for these parameters are discussed in Section 5.3 and compared with the results of the literature.

Observations by Sideband Digital Holography

Experimental setup

Holographic setup

The experimental setup is shown schematically in Fig. 1. A laser beam, with wavelength λ = 650 nm (angular frequency ω L ) is split into a local oscillator beam (optical field E LO ) and an illumination beam (E I ); their angular frequencies ω LO and ω I are tuned by using two acoustooptic modulators (Bragg cells with a selection of the first order diffraction beam) AOM1 and AOM2: ω LO = ω L + ω AOM1 and ω I = ω L + ω AOM2 , where ω AOM1,2 ≃ 2π × 80 MHz. The first beam (LO) is directed via a beam expander onto a CCD camera, while the second beam (I) is expanded over the surface of the reed, which vibrates at frequency f . The light reflected by the reed (field E) is directed toward the CCD camera in order to interfere with the LO beam (E LO ). 4 phases were used (phase shifting digital holography) and we select the first sideband of the vibrating reed reflected light by adjusting ω AOM1,2 to fulfil the condition:

ω AOM1 -ω AOM2 = 2π(f + f CCD )/4
, where f CCD is the CCD camera frame frequency. The complex hologram signal H provided by each pixel of the camera, which is proportional to the sideband frequency component of local complex field E, is obtained by 4-phases demodulation: H = (I 0 -I 2 )+j(I 1 -I 3 ) where I 0 . . . I 3 are 4 consecutive intensity images digitally recorded by the CCD camera, and j 2 = -1. From the complex hologram H, images of the reed vibration are reconstructed by a standard Fourier holographic reconstruction calculation [START_REF] Schnars | Direct recording of holograms by a CCD target and numerical reconstruction[END_REF]. These holographic reconstructed images exhibit bright and dark interference fringes. Counting these fringes provides the amplitude of vibration of the object (in the direction of the beam), which depends on the wavelength λ of the laser, and on the first Bessel function J 1 , for instance ±95nm for the first, ±770nm for the 5 th and ±1.6µm for the 10 th maximum (bright fringes) [START_REF] Joud | Imaging a vibrating object by Sideband Digital Holography[END_REF][START_REF] Joud | Fringe-free holographic measurements of large-amplitude vibrations[END_REF] 2 . This method has 3 main advantages: (i) The time for data acquisition is very short, about 3 minutes for recording 184 holograms, including holographic reconstruction.

(ii) The signal to noise ratio is significantly better than with traditional technology, particularly through the elimination of signal at zero frequency.

(iii) The visualization of large-amplitude vibration (order of magnitude: 0.1 mm) is possible by using high harmonics orders (up to several hundred times the excitation frequency).

Reed excitation

The reed was excited by a tweeter loudspeaker screwed onto an aluminium plate, connected to a clarinet mouthpiece. The lay of this mouthpiece was modified to be strictly flat. A plastic wedge of uniform thickness has been inserted between the lay and the reed, longitudinally to the same height as the ligature (Vandoren Optimum), allowing free vibrations of the entire vamp (length: about 38 mm), see Fig. 1. This ensures precise boundary conditions, avoiding any dependence to deformations of the reed. The repeatability of the longitudinal placing of the wedge and of the reed was ensured by a Claripatch ring [START_REF]Claripatch SA website[END_REF].

This setup requires some comments: (i) The reed is excited exclusively through the bore of the mouthpiece.

(ii) The pressure field in the chamber of the mouthpiece was not measured. Like for a real instrument, the edges of the reed (protected by the walls of the chamber) are subject to a pressure field, which is probably lower than the pressure acting on the rest of the vamp.

(iii) The boundary conditions are very different from those of a real instrument (no curved lay, no contact with the lip). In addition, the reed was not moistened for the measurement.

Figure 2: Typical holographic patterns of series A (asymmetrical sinusoidal excitation: the left side of the reed is more strongly excited than the right side). Frequency range: 1.4 to 20 kHz by steps of 25 cents (181 pictures ordered from left to right, continued on the next row; the tip of the reed is down on each picture). Some modes are easily identified: F 1 (1st row, 19th picture), T 1 (3rd row, 10th picture), F 2 (4th row, 7th picture), T 2 (4th row, penultimate picture), X1 (5th row, 5th picture) etc ... Modes T 3 and F 3 are almost at the same frequency (6th row, 4th and 5th pictures, probably). The last picture of the 3rd row corresponds to an acoustic resonance of the excitation device. It is present on all holograms of both series at the same frequency (examine Fig. 3). The excitation amplitude exponentially increases until the 73rd picture (3rd line, 21th picture), being hold constant beyond.

(iv) The excitation device is almost closed. The acoustical resonances of the excitation device are unknown, but may quite easily be deduced by comparing different measurements, because they are always present at the same frequency.

Observation of resonance frequencies

Experimental protocol

55 clarinet reeds of model Vandoren V12 were purchased in a music shop: 12, 12, 20 and 11 reeds of nominal forces 3, 3 1 2 , 4 and 4 1 2 , respectively. 29 reeds were used for two preliminary studies in order to develop the measurement protocol. Each of these reeds was played a total of some tens of minutes, spread over several weeks before measurement with the final protocol. The other 26 reeds were strictly new by measurement, which was performed immediately after package opening (for 21 of them with the new hermetically sealed package by Vandoren, ensuring a relative humidity between 45 and 70%, according to the manufacturer), without moistening Figure 3: Typical holographic patterns of series B (symmetrical sinusoidal excitation). The modes X1 and T 3 cannot be distinguished anymore. Notice that T 1 is less marked than under asymmetrical excitation (for some reeds even difficult to identify) and that the pattern has a significant flexural component, strongly dependent of the lateral placing of the reed on the mouthpiece. Notice that the symmetry of the patterns near T 1 depends on the excitation frequency. the reed. Each reed was subject to 2 series of measurements:

• Series A (asymmetrical excitation: see Fig. 2): the right half of the mouthpiece chamber was filled with modeling clay to ensure a good excitation of antisymmetrical modes. 184 holograms were made ranging from 1.4 to 20 kHz (sinusoidal signal), by steps of 25 cents. The amplitude of the excitation signal was exponentially increased in the range 1.4 to 4 kHz, from 0.5 to 16 V, then kept constant at 16 V up to 20 kHz. This crescendo limits the amplitude of vibration of the first two resonances of the reed. The temperature was not measured (about 20 • C).

• Series B (symmetrical excitation: see Fig. 3): the modeling clay was removed. The protocol is otherwise identical to this of the first series. The reeds were inadvertently exposed during one night to the very dry and warm air from the optical laboratory between the two series of measurements. The reeds lost between 2 and 4% of their mass. In what follows we try to interpret the influence of this fact. The temperature was around 23-25 • C.

Nomenclature of normal modes

Distinguishing 3 morphological classes, we classify the modes of a clarinet reed as follow : i) The "flexural" (or "bending", or "longitudinal") modes, listed below F , whose frequencies mainly depend on the longitudinal Young modulus (E L ) and polarized mainly in the z axis, ii) the "transversal" (or " torsional", or "twisting" ) modes, listed below T , mainly dependent on the shear modulus in the longitudinal / tangential plane (G LT ), and iii) the "generic" (or "mixed") modes, listed below X, sensitive to both moduli E L and G LT (see Fig. 4). A subclass of flexural modes may be distinguished: the "lateral" modes (listed below L), polarized mainly in the y axis (see Fig. 5). These modes were not observed in our study. In the software Catia used for the simulations, the orthotropic model is Cartesian and not cylindrical. Therefore an exact equivalence between x, y, z, and L, T, R respectively can be assumed since we observed no important deviation between the direction of the grain and the axe of symmetry of the reed.

The modes have been numbered after the order of increasing frequencies from a preliminary modal analysis we performed. In our analysis, however, the identification of a mode is based upon morphological criteria. As a matter of fact, the mode number and the order of observed frequencies are not necessarily identical for all reeds.

Strictly speaking, the optical method only allows to observe the resonance frequencies of the reed and not the eigenfrequencies. Therefore the observed deformation patterns are a priori not identical to the eigenmodes of the reed. Nevertheless in practice no major differences have be found between the computed eigenmodes (see Section 4) and the observed or computed deformation for a forced asymmetrical excitation at the corresponding frequency. For this reason we use the terminology "mode" for the maximum amplitude of the response of the reed to a forced excitation. This is somewhat abusive, because the small shift between the resonance frequencies due to damping and the eigenfrequencies computed by FEM, without damping, is ignored. Besides damping, the acoustic load is also able to shift the resonance frequencies. We assume that this discrepancy is approximatively the same for all reed.

Analysis of holograms; mode identification

More than 30000 holograms were made for this study and analyzed as follows: The picture where the number of interference fringes is locally maximum is determined. For some cases, we chose the hologram that is most similar to our numerical simulations (by Finite Element Method, see Section 4) or to other holograms (see Fig. 6). The holograms corresponding to an acoustical resonance of the system, present at the same frequency (4309 Hz) for all reeds, have been eliminated. The identification of the different patterns to those calculated by FEM was often quite simple. An exception have been encountered for F 3 and T 3, whose frequencies were often so close that our identification is sometimes uncertain. More sophisticated techniques would certainly solve this problem. Notice that other boundary conditions (e.g. with clamping closer to the tip of the reed) would also easily separate these two modes. The frequency of some higher modes could not always be measured, either because their frequency was beyond 20 kHz, or because their pattern could not be clearly identified.

Statistical analysis of resonance frequencies

4 flexural, 5 transversal, and 6 generic modes have been identified, namely all 15 first modes of the reed, excluding lateral modes. This number is significant, compared to the 4 modes detected by Pinard et al. [START_REF] Pinard | Musical quality assessment of clarinet reeds using optical holography[END_REF]. The 6th mode (L1) could not be identified, as it is a lateral mode (flexural mode moving mainly in the y axis), not excited by our loudspeaker. We tried to observe it by rotating the mouthpiece to the side, without success. Notice that higher modes could probably be identified using an ultrasonic loudspeaker.

Statistics

The statistics are displayed on Fig. 7 and detailed in Appendix A, with the analysis of correlations. For 14 measurements of resonance frequencies of the two series, identification of the total number of reeds (55) have been done. For other measurements, identification has been done only for a part of this number. The value of the ratio of the standard deviation σ to the mean value µ, i.e. the relative standard deviation, is found to be between 2 and 5% (about 1/3 tone). If we admit Gaussian distribution for the measured frequencies, 99% of the observations typically range about ±1 tone (±200 cents) around the mean value (i.e. µ ± 3σ), for all frequencies.

The identification of the mode X6 is uncertain: it seems to appear for frequencies lower than those of our simulations. Mode T 5 is on the limit of the range we studied: this explains the small value of the standard deviation.

Between series A and B, the flexural modes F 1 to F 4 lower their mid range, while the transversal modes slightly increase it. The difference between the two series probably lies mainly in the drying of the reeds, and this seems to have a statistically significant effect. This is surprising, because drying decreases the density of the reed, and theoretically this should proportionally increase all frequencies. In addition, according to Obataya et al. [START_REF] Obataya | Acoustic properties of a reed (Arundo donax L.) used for the vibrating plate of a clarinet[END_REF], drying is expected to increase E ′ L (at least around 400 Hz), which should also increase the resonance frequencies. However Chevaux [START_REF] Ph | Améliorations de la durée de vie des anches doubles en roseau pour instruments à vents[END_REF] observed that drying diminishes E ′ L for material extracted from the inner side of the cane and augments slightly E ′ L for material extracted nearer from the outer side (for cane suitable for oboe reeds), at least in the frequency range 100-500 Hz.

The hypothesis of an influence of the excitation method on the resonance frequencies seems unlikely, as well as the hypotheses of a poor reproducibility of the position of the reed on the mouthpiece between measurements or of the modification of the acoustic load, due to the modeling clay.

Principal component analysis

Principal Component Analysis (PCA) is mathematically defined as an orthogonal linear transformation transforming the data to a new coordinate system, such that the greatest variance by any projection of the data comes to lie on the first coordinate (called the first principal component or first factor), the second greatest variance on the second coordinate, and so on [27]. Theoretically PCA is the optimum linear transform for given data in terms of least squares. PCA is based upon the calculation of the eigenvalue decomposition of the covariance (or of the correlation) matrix (see e.g. [START_REF] Jolliffe | Principal component analysis[END_REF]).

A PCA has been performed using the FACTOR module of SYSTAT [28]. The 14 variables (observed frequencies) presenting complete measurements for all reeds have been selected (all variables having 55 identified pattern, N A or N B = 55, see Table 7). Frequencies are rated in cents.

The 4 largest eigenvalues have been selected. They capture 91.2% of the total variance of our sample (respectively 53.6%, 21.4%, 10.8% and 5.4% for each factor). A fifth factor would

AF1 BF1 AF2 BF2 AF3 BF3 AF4 BF4 AT1 BT1 AT2 BT2 AT3 BT3 AT4 BT4 BT5 AX1 BX1 AX2 BX2 AX3 BX3 AX6 BX6 BX4 BX5 G6 C7 E7 G7 C8 E8 G8 C9 E9 G9 C10 E10
Figure 7: Observed resonance frequencies for the different reeds of each series, according to the results of Table 7. AF1 signifies mode F1, series A. The series A is in black, and the series B in white (yellow online). In abscissa the frequencies in logarithmic (musical) scale. From the left to the right, the different notes correspond to the following frequencies: capture only 2.5% variance more. The 14-dimensional data have been linearly projected onto a 4-dimensional factor space.

The factor space can afterwards be orthogonally rotated, for instance for maximizing the correlations between rotated factors and observed variables. In the studied case, no a priori knowledge about the orientation of the factor space is available. For an easy comparison, using the VARIMAX algorithm, we choose to maximize the correlations between rotated factors and all available variables (observed resonance frequencies and theoretical components from the model described hereafter in Section 5.2).

We performed also a PCA separately for each measurement series (A and B: 9 and 5 variables, respectively). From series A we detected 3 important factors capturing 90.8% of the variance (56.9, 23.0, 11.0%, respectively). From series B we detected also 3 important factors capturing 94.1% of the variance (54.0, 26.8, 13.3%, respectively). A 4th factor would capture only 3.6% more for series A and 3.4% for series B. One factor seemingly disappeared, compared with the PCA performed on both series. A hypothesis is that this factor is related to the hygrometric change between the two series.

>From Table 1, we see that all transversal and generic modes are well correlated with f actor 1 ; f actor 2 correlates with high frequency flexural modes of both series (however notably better with those of series A); f actor 3 well correlates with F 1 of both series (and somewhat with other low frequency modes: AT1, BT1 and AF2), whereas f actor 4 correlates quite well with high frequency flexural modes of series B.

Conclusions from the statistical analysis

Fifteen modes of vibration of the clarinet reed have been observed, while previous studies investigated 4 to 5 modes only [START_REF] Guimezanes | Etude Expérimentale et Numérique de l'Anche de Clarinette[END_REF][START_REF] Pinard | Musical quality assessment of clarinet reeds using optical holography[END_REF]. The observed resonance frequencies are often highly correlated, especially those among the "transversal" modes and, to a lesser extent, those among the "flexural" modes. The nominal reed strength is surprisingly better correlated with the frequencies of "transversal" modes as with those of "flexural" modes. The flexural modes within the same series are poorly correlated.

A principal component analysis of the resonance frequencies identifies 4 main factors, capturing 91.2% of the variance of the sample. The data can therefore be reconstructed with 4 uncorrelated factors only (error: RM SD = 21.8 cents, see Appendix D.2 and E). The effect of hygrometric change between both measurement series can seemingly be described with 1 factor only.

These statistical facts offer a guidance for modeling appropriately the mechanics of the clarinet reed.

Development of a mechanical model 4.1 Choice of a viscoelastic model

In the present study our concern is to develop a model with a minimal number of physically related components, that adequately reconstructs the observed resonance frequencies of our reeds. We presume that these components offer an objective characterization of the material composing each reed. A sensitivity analysis by FEM calculation assuming an orthotropic, elastic material has been conducted, and showed that the longitudinal Young modulus E L and the longitudinal / transverse shear modulus G LT play a leading role. Nevertheless taking into account the previous result of 4 factors given by PCA, we will see that an elastic model is not sufficient to establish a satisfactory model with 2 degrees of freedom only (i.e. variables E L and G LT per series. Therefore a viscoelastic model is sought.

It is well known that the stiffness of natural materials like wood or cane varies with the frequency of the applied stress and with the temperature. The material is stiffer at low temperature and at high frequency. At low frequency or high temperature the material is almost perfectly elastic and reaches the rubbery modulus. At high frequency or at low temperature the glassy modulus is reached; the material is almost perfectly elastic, also, but stiffer. At mid frequency or mid temperature, the apparent modulus (called storage modulus, i.e. the real part of the complex Young modulus for this frequency) is between the two values. For a particular frequency, called relaxation frequency, the storage modulus is exactly at the average of glassy and rubbery moduli. Around this frequency dissipation is maximum. Once the characteristic curve is known (for given temperature and different frequencies, or for given frequency and different temperatures), the Arrhenius equation 3 offers usually an adequate estimate of the stiffness for any frequency and any temperature, within a quite broad range [START_REF] Roylance | Engineering viscoelasticity in 3.11 OpenCourseWare[END_REF].

The determination of the mechanical parameters of a natural material requires to determine for each axis of the orthotropic material the value of 3 parameters (Young modulus, shear modulus and Poisson's ratio). These 9 parameters may exhibit viscoelastic behavior, requiring theoretically for each one the fit of a viscoelastic model, such as the general linear solid (also called Zener model or 3-parameter model, see [START_REF] Roylance | Engineering viscoelasticity in 3.11 OpenCourseWare[END_REF][START_REF] Chaigne | Acoustique des instruments de musique[END_REF][START_REF] Dill | Continuum mechanics: elasticity, plasticity, viscoelasticity[END_REF][START_REF] Gaul | Experimental Determination and Modeling of Material Damping[END_REF]) 4 . The chosen model is based on 3 The shift in relaxation time is:

Ln(shif t) = Ea R 1 T -1 T ref
, where Ea is the activation energy, R is the gas constant (8.314 J/K mol) and T and T ref the absolute temperatures in K. For instance, a shift of +10 • C from a reference temperature of 20 • C decreases the relaxation time by 16%, if Ea is 13 kJ/mol. 4 Other multidimensional viscoelastic models could be also considered. In order to fit a wide range of frequencies (more than 2 decades), a 4-parameter model with fractional derivative would be required [START_REF] Gaul | Experimental Determination and Modeling of Material Damping[END_REF]. In addition, these parameters are known to be sensitive to moisture content. Moreover, the cane is not homogeneous. The stiffness varies in radial direction [START_REF] Ph | Améliorations de la durée de vie des anches doubles en roseau pour instruments à vents[END_REF] and local irregularities may be important, as shown by J.-M. Heinrich [START_REF] Heinrich | Recherches sur les propriétés densitométriques du matériau cane de Provence et ses similaires étrangers; relation avec la qualité musicale; étude associée d'une mesure de dureté[END_REF].

6 parameters (3 parameters for both variables E L and G LT ). Therefore the viscoelastic model (Section 4.4) has many degrees of freedom (6, for each of the 2 series of measurements, compared to the 4 factors detected by PCA for both series), for solving adequately the inverse problem (see Section 5 for the reduction of the number from 12 to 4).

Computation method

Considering viscoelasticity leads to complex modes with complex eigenfrequencies. Compared to the non-dissipative, elastic case computed by FEM, the main consequence of viscoelasticity, besides dissipation, is that stress and strain are not in phase. For sake of simplicity, we limit the computation to eigenfrequencies only, and assume that they depend on the storage moduli only (i.e. dissipation has a negligible influence). Having reduced the viscoelastic problem to an associated elastic one, the elastic solution may be used (see e.g. Ref. [START_REF] Roylance | Engineering viscoelasticity in 3.11 OpenCourseWare[END_REF]). In order to compute the resonance frequency ω r after an elastic model, according to Ref. [START_REF] Gaul | Experimental Determination and Modeling of Material Damping[END_REF], we admit that E ≃ E ′ (ω r ), where E ′ (ω r ) is the real part of the complex modulus in the frequency domain. This hypothesis implies that the calculation of the eigenfrequencies from the values of the storage modulus is done by an iteration procedure, and allows to use a FEM software (Catia) which does not allow computing with frequency-dependent coefficients.

Therefore we first present results of FEM simulations (Section 4.3), assuming an elastic and orthotropic behavior of the reed (modeled in Section 4.3.1). This helps to identify the modes in experiments, and allows obtaining a fit formula (Section 4.3.3) for computing the 11 lower resonance frequencies with respect to two parameters only, E L and G LT , detected after a sensitivity analysis (Section 4.3.2). The fit formula, called "metamodel", is then used in the iterative procedure for the computation of the viscoelastic model. It allows a great reduction of computation time, compared to the FEM, and this is very useful for the inverse problem. Such a metamodel could be directly computed for the viscoelastic model with an appropriate software, but starting with the elastic model simplifies the fitting procedure.

Elastic model

Modeling the reed

The clarinet reed is defined in a Cartesian axis system x, y, z (see Fig. 5). The origin is located in the bottom plane, at the tip of the reed. The material is defined as 3D orthotropic and assumed to be homogeneous, whose longitudinal direction L is parallel to the x axis, the tangential direction T parallel to the y axis and the radial direction R parallel to the z axis 5 .

The dimensions in the xy plane are consistent with the measurements given by Facchinetti et al. [START_REF] Facchinetti | Numerical and experimental modal analysis of the reed and pipe of a clarinet[END_REF]. The heel of the reed is made out of a cylinder section, diameter 34.8 mm, maximum thickness 3.3 mm. The shape of the reed is defined in Appendix B.

During playing, the reed has two contact surfaces with the ligature. For the present simulations, the reed is clamped in the same way than for normal playing, on two rectangular surfaces 23 × 1 mm, spaced laterally by 5 mm, 38.2 mm from the tip of the reed, simulating the contact surfaces on the Vandoren Optimum ligature. However, unlike normal playing, the whole vamp of the reed is free to vibrate (see Figure 1).

For the simulations, the "Generative Part Structural Analysis" module by Catia v.5.17 (Dassault Technologies) is used, with mesh Octree3D, size 2 mm, absolute sag 0.1 mm, parabolic tetrahedrons. The generated mesh involves 5927 points, allowing both a good accuracy and a reasonable computing time (around 35 seconds).

Coefficient

Sensitivity analysis of elastic coefficients

For selecting the most relevant parameters, we conducted a One-At-a-Time sensitivity analysis [START_REF] Saltelli | Global Sensitivity Analysis[END_REF], varying each coefficient by ±10% and computing the first 16 modes, based on the following reference values:

E L =14000 MPa, E T = E R = 480 MPa, ν LT = ν LR = ν T R =0.22, G LT = 1100 MPa, G LR = G T R = 1200
MPa. The density ρ was set to 520 kg/m 3 , according to the estimation by Guimezanes [START_REF] Guimezanes | Etude Expérimentale et Numérique de l'Anche de Clarinette[END_REF]. The results are shown in Table 2. Notice that E L and G LT plays a decisive role, while E T plays a marginal role and all other parameters have an almost negligible influence on the resonance frequencies. As a consequence, the moduli E L and G LT are the variables retained in the model. The approximate value of E T has been estimated according to the morphology of the patterns of higher order modes. This value is consistent with measurements given by Spatz et al. [START_REF] Spatz | Biomechanics of the giant reed Arundo donax[END_REF].

Notice that these results show the validity of a 2D approach, the reed being modeled as a thin plate. This should be used for further studies.

Metamodel approximating the resonance frequencies

The following analytic formula ("metamodel") predicts quickly and efficiently the resonance frequencies of a clamped/free clarinet reed. It was established in the following way: Frequencies of the first 16 modes were computed by FEM, according to a network of 92 separate pairs of values for E L and G LT , ranging from 8000 to 17000 MPa and 800 to 1700 MPa, respectively. The other elastic coefficients were held constant, according to the reference values cited above. For the range of simulation values, this arbitrary formula (developed by trial and error) provides a very good fit (generally better than ±5 cents, see Table 9). Expected resonance frequencies f are first found in cents (F C) from the note F6 (1396.9 Hz), and finally in Hz: The index m is the number of the mode defined in Appendix C, Table 9, where the values of the coefficients a m,q are given (E L and G LT are expressed in MPa).

f (m, E L , G LT ) = 1396.9 × 2 F C/1200 , where (1)
F C = a m,0 + a m,1 E p + a m,2 G p + a m,3 E p G p + a m,4 E 2 p + a m,5 G 2 p , E p = E L -0.66643 and G p = G LT 0.7627 . E 2 E 1 E 3
The influence of the density is easy to predict: frequencies vary proportionally to ρ -1/2 . The computing cost of this metamodel is about 10 7 times lower than with FEM, largely simplifying the inverse problem.

Efficiency of the metamodel

Equation ( 1) can be used to estimate the values of E L and G LT , providing a faithful reconstruction of the observed resonance frequencies. Theoretically these values could be computed for any pair of modes, after their respective observed frequencies. Unfortunately, this method gives no consistent results. A least squares fit is a more robust technique for such a computation. This leads however to systematic errors in the predicted frequencies: low-order modes are systematically overestimated, while high-order modes are underestimated. This can be corrected by adjusting the coefficients a m,0 (from Table 9), but this cannot explain the bad correlation among flexural modes within the same series (see Table 6). According to the elastic model, these correlations should be in all cases greater than 0.998. A hypothesis for resolving this contradiction is that the moduli are varying with the frequency in an individual way for each reed. Thus in the next paragraph we consider a viscoelastic model, where E L and G LT are frequency dependent. This leads to the addition of some parameters, which are to our mind more important that the other elastic coefficients. The fit of such a model requires many observations at different frequencies, in order to reduce the influence of measurements errors and of local irregularities in the structure of cane.

Alternative hypotheses could be considered in this context, as damping, acoustic load [START_REF] Facchinetti | Numerical and experimental modal analysis of the reed and pipe of a clarinet[END_REF], local variations in stiffness or in density, local deviations in thickness, compared to the assumed theoretical model. However, these hypotheses are probably unable to explain the hygrometricinduced individual variations we observed for each reed, thus our preference for the viscoelastic hypothesis.

Viscoelastic model

In this section, a Zener model is considered (see e.g. Refs [START_REF] Roylance | Engineering viscoelasticity in 3.11 OpenCourseWare[END_REF][START_REF] Chaigne | Acoustique des instruments de musique[END_REF][START_REF] Dill | Continuum mechanics: elasticity, plasticity, viscoelasticity[END_REF][START_REF] Gaul | Experimental Determination and Modeling of Material Damping[END_REF]). This model is applied to both moduli E L and G LT . The scheme of the standard viscoelastic solid is presented on Fig. 8, with two springs E 1 and E 2 and a dashpot E 3 6 . At low frequencies, E 2 and E 3 have practically no effect (the rubbery modulus E 1 dominates). At high frequencies, E 3 has practically no effect (the glassy modulus E 1 + E 2 dominates). In the frequency range near E 2 /(2πE 3 ) the dissipation due to E 3 is maximal and the apparent modulus (storage modulus) is in the mid-range. The stress σ and the strain ε are related by the constitutive equation:

σ + τ 1 σ = E 1 (ε + τ 2 ε) (2)
in which τ 1 = E 3 /E 2 is called the relaxation time and τ 2 = E 3 (E 1 + E 2 )/(E 1 E 2 ) the retardation time. E 1 is called rubbery modulus and E 1 + E 2 glassy modulus. In harmonic regime, for an angular frequency ω, the Young modulus is complex:

E * (ω) = E 1 + E 2 - E 2 2 E 2 + jωE 3 = E 1 + E 2 + -E 2 + jωE 3 1 + (ωE 3 /E 2 ) 2 (3)
The second formulation separates the real part (E ′ (ω): storage modulus) and the imaginary part (E ′′ (ω): loss modulus) of E * (ω). The storage modulus can thus be written as:

E ′ (ω) = E 1 + E 2 - E 3 2 E 2 2 + ω 2 E 2 3 = E 1 1 + ω 2 τ 1 τ 2 1 + ω 2 τ 2 1 . ( 4
)
Notice the properties:

E ′ (0) = E 1 ; E ′ (1/τ 1 ) = E 1 + E 2 /2 E ′ (∞) = E 1 + E 2 ; ∂E ′ ∂ω 1 τ 1 = E 3 2 .
For the sake of simplicity, the parameters corresponding to E L and G LT are denoted E 1 , E 2 , E 3 , and G 1 , G 2 , G 3 , and the storage moduli given by Equation (4) E ′ (ω) and G ′ (ω), respectively. Therefore for each reed (and each series), the model requires 6 parameters instead of 2 (while experiments gave 4 main factors only for the whole set of results).

>From the knowledge of the 6 parameters, the resonance frequencies can be deduced by an iteration procedure. For each mode the starting point of the iteration is the mean value f (0) of the experimental resonance frequency (see Table 7), then the storage moduli are deduced from Equation (4), then a new value f (1) by using Equation (1), etc... The convergence of the iteration method is fast, actually one iteration is enough. This can be understood by the fact that the derivative of the iterated function is small (notice that the two first rows of Table 2 correspond to the derivative of E L and G LT with respect to frequency). If we give an arbitrary value, for instance f (0) = 6000, one iteration more is required for a comparable precision. In all hypotheses (see Appendix D.6), we used one iteration only. This procedure allows the determination of f (m) from the coefficients E 1 , E 2 , E 3 , G 1 , G 2 , G 3 for a given reed and a given series.

Inverse problem and selection of a robust model

Simplification of the model by multiple regression

In order to solve the inverse problem for each reed (and each series), we use a classical Mean Squared Deviation method, from the experimental values of the 11 resonance frequencies listed in Table 9 (see Appendix C). The results are given and discussed in Section 5.3.

However for each reed the viscoelastic model provides 12 parameters, i.e., 12 DoF, and the PCA showed that this number needs to be reduced. Actually this model conducts sometimes to non-physical results (negative rubbery modulus, for instance). This problem comes out because the observed resonance frequencies are far from 0, so the rubbery modulus E 1 cannot be estimated precisely. For that purpose multiple regression (see Appendix D for details) is used, by introducing correlations among parameters and reducing the degrees of freedom to a number of 4, called "components", which are linearly related to the parameters. The 4 components are very similar to the 4 factors computed by the PCA, but, because Eqs. (1 and 4) are nonlinear, a small deviation is inevitable for optimal results. Factors and components are consequently strongly correlated (>0.95, see Table 1).

We tested different hypotheses to establish a satisfactory robust model (denoted H1 to H9 and described in Appendix D), together with the detailed computation method. Each hypothesis leads to a given number of components related to the parameters through the regression coefficients. Notice that the parameters and components depend on both the reed and series, while the regression coefficients which correlate the parameters are independent of the reed and the series. A constant value elastic model (the parameters are fixed, and do not depend on the reed or the series, Hypothesis H1) is not sufficient. Similarly a 2-parameter elastic model (E L and G LT , independent of the series, Hypothesis H2) is not sufficient as well. Conversely a linear model without constraints (with 60 independent coefficients, corresponding to a matrix of order 12(4+1) elements7 , Hypothesis H7) is not necessary, because many coefficients are very small. Eventually a 9 coefficient linear model (H4) has been found to be very satisfactory and is described in Section 5.2. The following ideas were applied: after eliminating the very small coefficients and fitting the observations with the remaining coefficients, it is observed that the predictive quality of the model is almost not affected by the simplification. This was done step by step. The problem of non-physical values for certain parameters can be corrected by setting that the damping parameters E 3 and G 3 are constant, independent of the reeds and of the series (typical values: E 3 = 0.28 and G 3 = 0.02). It remains 4 parameters for each series and reed, E 1 , E 2 and G 1 , G 2 (i.e. 8 parameters for each reed). This ensures generally that these parameters fall in a plausible range, when fitting the model. Moreover a hierarchical structure can be introduced in the model, isolating the hygrometric component, bringing the remaining 3 components to a common basis (Section 5.2) and simplifying the problem (reduced to only 9 regression coefficients) and giving some insight in the data structure.

The RM SD (Root Mean Square Deviation, see Appendix D.2) is found to be 30.4 cents for Hypothesis H4, very close to 29.8 cents for H5 with 9 coefficients more. Moreover the standard deviation of the residuals for Hypothesis H4 (and also Hypothesis H5) varies very few over the different resonance frequencies (all around 30 cents).

The fit quality cannot be considered as a perfect and definitive proof that our model reflects the true values of the corresponding storage moduli. The influence of some missing parameters in the model should be examined (for instance differences in thickness between reeds, non constant modulus E T , non constant density ρ or radial variation of E L ). Anyways, the presented model reflects real mechanical differences between the reeds, very similar to those objectively detected by the PCA.

Robust estimation of the parameters of the viscoelastic model (Hypothesis H4)

Hypothesis H4 is chosen so that no coefficient can be removed without impacting notably the quality of fit. It can be thought as the minimal structure allowing an adequate reconstruction of the observed resonance frequencies, in conjunction with the viscoelastic model (Eq. ( 4)) and the metamodel (Eq. ( 1)). This minimal structure makes the model more robust against measurements errors, even if it probably introduces some bias. As a first step, our concern is to eliminate the influence of the moisture content and to bring both series of measurements to a common basis (i.e. predict the effect of drying on the viscoelastic parameters of series B, the series A being taken as a reference). For series A, the components remain unmodified (series A is taken as reference):

ě1 [s = 1, n] = e 1 [n] ě2 [s = 1, n] = e 2 [n] ě3 [s = 1, n] = e 3 [n] (5)
For series B, the effect of drying on the components is predicted by:

ě1 [s = 2, n] = c 10 + e 1 [n] ě2 [s = 2, n] = c 20 + 1 2 (e 2 [n] + e 4 [n]) ě3 [s = 2, n] = c 30 + e 3 [n] (6)
With this choice of components, the viscoelastic parameters of the model for series s and reed n can then be estimated as follows:

E 1 = d 10 + d 12 (ě 2 [s, n] + ě3 [s, n]) E 2 = d 20 + d 12 (2ě 2 [s, n] -ě3 [s, n]) E 3 = d 30 G 1 = d 41 (6 + ě1 [s, n]) G 2 = 2d 41 (3 + ě1 [s, n] -ě3 [s, n]) G 3 = d 60 (7)
This implies some other interesting relationships: The change in density between the two series of measurements was not measured precisely (about -2 to -4%). In the model, density is considered as constant.

E 1 + E 2 = d 10 + d 20 + 3d 12 ě2 [s, n] G 1 + G 2 = d 41 (12 + 3ě 1 [s, n] -2ě 3 [s, n]) 2E 1 -E 2 = 2d 10 -d 20 + 3d 12 ě3 [s, n] 2G 1 -G 2 = 2d 41 (3 -ě3 [s, n]) (8)
Fig. 9 shows an approximation of the frequency dependence of E L and G LT , computed for the mean value of all reeds. For Series A, the storage modulus E ′ (ω) increases from 11700 MPa at 2 kHz (F 1) to 17200 MPa at 16.8 kHz (F 4), while for Series B, it increases from 11100 to 13300 MPa. Therefore under "normal" hygrometry the reed bends in a notably viscoelastic manner, whereas the ultra-dry reed bends in an more elastic manner. For G LT , the reed is generally notably viscoelastic, according to our model. The corresponding values are: 783 to 1290 MPa for Series A, and 896 to 1436 MPa for Series B. Corresponding statistics are displayed on Table 3. Drying seems to increase G LT and decrease E L (except around 2300 Hz), explaining the good correlation between the variables AF1 and BF1.

This simplified model permits interesting conclusions about the structure of our data:

• The component e 1 is related exclusively to G LT ;

• e 2 is related exclusively to E L ;

• e 3 increases proportionally to the rubbery modulus of E L (and not its glassy modulus), and decreases proportionally to the glassy modulus of G LT (and not its rubbery modulus);

Increase in

E L G LT storage
-E ′ (2πf F 1 )/E ′ (2πf F 4 ) or 1 -G ′ (2πf F 1 )/G ′ (2πf F 4 ), with f F 1 = 1996Hz
and f F 4 = 16784Hz. Between both series, the correlation is 0.83 for E L and 1 for G LT (as a consequence of the simplification of the model) .

it regulates therefore the viscous component common to E L and G LT8 ;

• e 4 takes into account the variation in moisture content between series A and B.

Notice that the product of the rubbery moduli for the series A, E 1 G 1 is correlated at 72% with the nominal reed strength. Individually, these moduli correlate only at 50% and 55% with the nominal reed strength, respectively.

Results and discussion

Tables 4 and5 show our results for hypotheses H4 and H5. A comparison with results by other authors is difficult or even quite impossible, because of the disparate structure of the measurements. Such a comparison requires the reconstruction of the measurement data (when possible), the fit of a viscoelastic model or the extrapolation of the values, in order to reach the frequency and temperature range of our measurements. The validity of such a highly speculative task is questionable. For the storage modulus E ′ L , all reconstructed values from other authors fall in the range around the average of our measurements ±3 times the standard deviation, however most of the time in the lower range. This probably shows that the selected value for the density ρ is somewhat too high. No representative statistics are available by the other authors.

The most important disagreement, compared with our model (hypotheses H4 and H5), is about the relaxation frequency. The explanation is probably because the studied frequency range was not same. For G LT , we found no viscoelastic measurements by other authors. Our results are summarized in Tables 4 and5. For E ′ L , between hypotheses H4 and H5, f r and E ′ L at 4 kHz agree well, whereas E 1 and E 2 diverge by about 1-2 SD. This divergence comes because the observed frequency range was not broad enough. For the shear modulus, G 1 and G 2 are in good agreement for both hypotheses (and consequently f r and G ′ L at 4 kHz also). Our model is valid only for "ambient dry" reeds (since the ultra-dry conditioning was not controlled), in a frequency range which should not exceed one decade. We checked that a fractional derivative model after Gaul et al. [START_REF] Gaul | Experimental Determination and Modeling of Material Damping[END_REF] is not necessary in our narrow frequency range. Such models are however really efficient to cover a broad frequency range. For instance, the data by Lord for dry material [START_REF] Lord | Viscoelasticity of the giant reed material Arundo donax[END_REF] could be fitted very well (E 0 = 8108 MPa, E 1 = 2964 MPa, p = 0.298 MPa•s α , α = 0.546). Notice that the order of the derivative (α) is 1 in our viscoelastic model. Our viscoelastic model is able to partially explain the bad correlations observed between flexural modes. In Table 6, the correlations are compared among observed resonance frequencies and computed modal frequencies, according to the viscoelastic model (hypothesis H4 and H5). It seems that additional hypotheses (such as an irregular thickness) should perhaps be considered for improving the model. However, we should remember that the determination of resonance frequencies are attached with some uncertainties, especially for the modes F 3 and T 3.

Model

Storage Young modulus E ′ L (ω) at about 20 • C E 1 [MPa] E 2 [MPa] f r [Hz] E ′ L at
In order to clarify this issue, let us examine if the residuals (observed resonance frequencies minus computed modal frequencies with hypothesis H59 ) contain some pertinent information. A PCA shows that perhaps 2 residual factors contain some interesting information (explaining 30% and 18% of the residual variance). The first residual factor is correlated with AT3 (0.86), BT2 (0.74), AT2 (0.73) and AX2 (0.66). All these modes depend strongly on G LT . An adjustment of the coefficient a m,2 from the metamodel for the transverse modes could probably cancel this systematic bias (remember that the coefficients a m,q are computed from a theoretical model, which is probably also biased). Indeed, an increase by 14, 21, 16 and 11% of this coefficient affecting the modes T 1, T 2, T 3 and T 4 makes the RM SD drop from 29.8 to 28.5 cents. The second residual factor is correlated with AF1 (0.36), AF2 (0.29), AT2 (-0.29), BF1 (0.28) and BT1 (-0.23). This probably reveals a competition between flexural and transversal modes when fitting the model. The bias probably comes from the coefficient a m,1 , regulating the linear dependance to E L in the metamodel. Adjusting the coefficients a m,1 for all modes and the coefficients a m,2 for all transverse and generic modes makes the RM SD drop down to 26.2 cents, reaching the size of the measurement steps (25 cents). The adjustments of a m,1 are very small for the flexural modes F 1 to F 4: -6, 1, 6 and 1%.

This shows that the most important bias depends linearly on the 2 most important parameters (E L and G LT ) of our FEM computations. No supplementary parameter is required until this bias is removed (theoretically down to a RM SD of 21.8 cents, according to H10). For hypothesis H4, the same linear adjustment of the metamodel lets the RM SD drop from 30.4 down to 27.0 cents.

Conclusion

The numerical model is satisfactory. From the statistical analysis are discussed in Section 3.3, it allows selecting the most important parameters describing the mechanical behavior of a reed.

The efficient elastic metamodel can be extended to a viscoelastic behavior of the reeds, approximating the resonance frequencies from the longitudinal Young modulus E L and the longitudinal / transverse shear modulus G LT and considering the hypothesis of their frequency dependence. A reconstruction of the observed resonance frequencies can be achieved with a good accuracy, estimating for each reed only 4 components, from which the parameters of a viscoelastic model are computed as a linear combination. The selected model (according to hypothesis H4) is probably slightly biased, but it is more robust against measurement errors than more refined models.

Table 1 shows that these components are highly correlated to the factors computed by PCA (0.96 to 0.98).

The proposed method allows the determination of 3 mechanical parameters characterizing the material composing each reed, with a single series of measurements, using Equations (1, 4 and 7). The reed should be conditioned with a relative humidity corresponding to the one ensured by the hermetically sealed package by Vandoren (about 55%). The fourth parameter cannot be determined in a reproducible way, since the exposure of the reeds to the ultra-dry air of the optical laboratory was not controlled. The same protocol and the same viscoelastic model can be used for other kinds of single reeds (bass clarinet, saxophone). Only the coefficients of Table 9 have to be recomputed after a FEM simulation of the corresponding reed shape.

Despite the fact that the eigenmodes of higher order probably play no important role in the acoustics of the clarinet, the present study shows that they reveal the inner structure of the material building the tip of the reed, so a new step could be done for an objective mechanical characterization of the clarinet reed. A subsequent study should examine if the obtained components are correlated with some musical qualities of the reeds. This could help the reed makers to gain a better control on their products. The thickness of the vamp at point (x, y) is interpolated as follows: first, we interpolate 3 points at y=0, 4 and 6mm, with 3 cubic splines, according to Table 8. These three points (s 0 (x), s 1 (x) and s 2 (x)) define a biquadratic polynomial: vamp(x, y) = p 0 (x) + p 1 (x) y 2 + p 2 (x) y 4 with p 0 (x) = s 0 (x), p 1 (x) = (-65 s 0 (x)+81 s 1 (x)-16 s 2 (x))/720 and p 2 (x) = (5 s 0 (x)-9 s 1 (x)+ 4 s 2 (x))/2880, allowing an interpolation on the y axis. The network of points above was estimated using a least squares fit, based on a network of 12 × 24 thickness measurements, achieved with a dial indicator and a coordinates-measuring table (estimated accuracy: ±5µm in z, ±50µm in x and y). We measured twenty reeds and select a particularly symmetrical one as reference. This method allows the reconstruction of the measurement network with an accuracy of ±10µm.

B Defining the shape of the reed

The thickness of the heel is defined to be :

heel(y) = -14.1 + 17.4 2 -y 2
The contour of the reed in the xy-plane is defined by:

contour(x) =          0 x <

D Development and selection of a simplified viscoelastic model

In this section we describe how the proposed model was developed and selected. Alternate options are presented.

D.1 Data structure

We need a specific notation for denoting our complicated multivariate data structure as arrays, after a list of indices (see Table 10). The 11 modes are defined after Table 9 (from 1 to 11: F 1, T 1, F 2, T 2, X1, T 3, F 3, X2, X3, T 4, F 4). We define a variable v n,s,i,j holding all parameters of our viscoelastic model: v n,s,1,j = E j and v n,s,2,j = G j for reed n and series s. (

The array r n,s,m holds the reconstructed resonance frequencies computed for all reed, series and modes, with the parameters v n,s,i,j of our viscoelastic model (see Section 4.4).

D.2 Mean Squared Deviation

As a cost function to minimize, we define the Mean Squared Deviation MSD (also called Mean Squared Error) between reconstructed and measured resonance frequencies o n,s,m :

M SD = 1 N SM N n=1 S s=1 M m=1 (o nsm -r nsm ) 2 (10)
With Eq. ( 10), the components of array v can be fitted by any appropriate algorithm for minimizing a multivariate function. All available measurements can be utilized for fitting the models 11 . Missing observations o m,s,n are then eliminated while computing M SD.

Our model allows a quite good reconstruction of the resonance frequencies, with a √ M SD ≡ RM SD (Root Mean Square Deviation) smaller than 20 cents (it is very small for lower modes, and always smaller than 25 cents). This corresponds to the hypothesis named H9 (all hypotheses H1 to H11 are presented and commented in Section D.6), but, as explained in section 5, the values of the coefficients are not always plausible physically. In Section D.3, we examine how to regulate this drawback by multiple regression.

D.3 Estimating the parameters of the viscoelastic model by multiple regression

Multivariate linear regression consists in projecting linearly a n dimensional space on a onedimensional space. The generic equation is

y = a 0 + N n=1 a n x n = N n=0 a n x n , if x 0 = 1 (11)
Equation ( 11) can be generalized for multiple regression as:

y = Ax (12)
where x and y are column vectors and A a matrix. In what follows, in order to use conventional vectors and matrices (with respectively one and two dimensions), we use the following notation:

v j [n, s, i] = v n,s,i,j is the jth element of the column vector v[n, s, i].
In our case, we have no prior knowledge about the relationships among the parameters of our model. We assume that each parameter in the model can be computed as a linear combination of some unknown independent components by multiple regression. The multiple regression formula can be written as:

v j [n, s, i] = K k=0 M jk [s, i] • e k [n] (13)
In conventional vectors and matrices notation, Eq. ( 13) reads:

v[n, s, i] = M [s, i] • e[n]
where e[n] is the vector of the orthogonal components for each observed reed n, M [s, i] is the regression matrix, independent of the reed number, depending on the series and the kind of modulus (E L or G LT ). v[n, s, i] is the vector of the parameters of the viscoelastic model.

We have only to choose an arbitrary number of components, for instance K = 4, referring to our PCA. We introduce arbitrary constraints in order to obtain more comparable results: over the different reeds, we state that the components must be orthogonally normalized (mean 0, standard deviation 1 and intercorrelation 0). The matrix of components e has consequently to satisfy:

e • e T =       N 0 0 0 0 0 N -1 0 0 0 0 0 N -1 0 0 0 0 0 N -1 0 0 0 0 0 N -1       (14)
Each individual column vector e[n] from this matrix is written as follows (for the component e 0 [n] = 1 : see Equation ( 11)):

e[n] = 1 e 1 [n] e 2 [n] e 3 [n] e 4 [n] T . ( 15
)
The components of M are fitted by minimizing M SD in Eq. [START_REF] Obataya | Mechanical relaxation processes due to sugars in cane (Arundo donax L.)[END_REF]. As starting value for the orthogonal components we set : e k [n] = f actor k [n], where f actor k [n] are the factors computed by PCA (see Section3.2). After a first estimation of M , it is possible to release the approximation about e: all components and all coefficients in the matrices can be fitted by the fitting procedure. However the number of variables to fit is probably much higher than allowed by the most algorithms of function minimizing. The fitting procedure has to be carried "by hand" with subsets of variables. The procedure we used is described below.

For K = 4 (the dimension of vector e[n] being 5), and some supplementary choices, the model works very well.

D.4 Empirical simplification of a 4-parameter model

We observed that the effect of hygrometric changes between both measurements series can be taken into account with one parameter only, practically without drop in quality of fit. This effect can be isolated on component e 4 [n] and the remaining components can be transformed linearly, so the further computations can be achieved from a common basis. Eq. ( 13) is then structured as:

M [s, i] = M [i] • M [s] , (16)
The reasoning considers the situation where the two series of observations are independent, with a number of components in vector ě

[s, n] = M [s] • e[n] reduced to Ǩ = K -1 = 3.
As a consequence the number of rows of matrix M [s] is 4, as well as the number of columns of M [i]. This approach, which allows separating hygrometry effects, offers a comfortable way to test different hypotheses, without changing the structure of the computation, by setting some coefficients in the matrices at some arbitrary values or by introducing some linear relationship between coefficients. We get fewer "active" coefficients to fit in the model: the fitting procedure is much faster and this raises the probability to find the best possible fit.

We tried to minimize the number of coefficients different from zero in the matrices, without substantial drop in quality of fit. The model was fitted "by hand", using the solver of Excel (Microsoft Office). We found empirically that a quite sparse setup still gives a good fit

This corresponds to the hypothesis named H5. Furthermore some coefficients may be proportional to others, without noticeable drop in quality of fit. This diminishes the number of active coefficients in the different matrices from 18 to 9 (hypothesis H4):

D.5 Adjusting coefficients for removing systematic errors

After fitting the different models, we observed some systematic deviations in the resonance frequencies between model and observations. This error has probably two different origins: an inevitable inaccuracy in the FEM computation (and in our metamodel) and an error for parameters not included in the model. A straightforward way to minimize the residuals is to fit the coefficients a m,0 in Eq. ( 1) 13 . Fitting all coefficients in Eq. ( 1) is doubtless a more questionable way to reduce this error (hypotheses H8 and H9). This can reduce the mean deviation between model and observations, but greatly increase the number of coefficients in the model (see however the discussion in Section5.3).

D.6 Testing different hypotheses

We tested different hypotheses with our model, in order to select a particularly efficient model. Some of them are summarized in the Table 11.

Hypotheses H1 to H3 (elastic model) present a poor fit; the adjustments for coefficients a m,0 are large, compensating partially for the missing viscoelastic components. All viscoelastic models are notably better and exhibit smaller adjustments for a m,0 . Frequency dependence for E L and G LT seems evident. Hypothesis H4 shows a good accuracy, with only 9 fitted coefficients (and 11 adjustments). Increasing the number of coefficients up to 60 brings only a marginal contribution (H5 to H7). Adjusting the other coefficients of Equation ( 1) in H8 and H9 improve the model, especially for the higher modes (Notice that no multiple regression is used for H9). Our FEM computations (and consequently our metamodel) are probably attached with systematic errors in this frequency range. The influence of E T should possibly be considered. Between H8 and H9, the total number of components (N × K) increases from 220 to 660. The adjustments within morphological classes are related: "flexural" modes shows systematically lower values than neighboring "transversal" modes.

D.7 Backward random validation

Following a suggestion by a referee, we applied a principal component analysis to simulated data computed after hypothesis H5: we assigned randomly a value for the 4 components and 55 reeds, following a normal distribution. We repeated this operation ten times. As expected, the PCA detected 4 factors capturing 91.2% (Standard Deviation 1.3%) of the variance of the simulated data (mean: 42.6, 20.9, 16.2 and 11.5% for each factor). This seems compatible with the analysis performed on the observed frequencies (Section 3.

E Reconstructing observed resonance frequencies by multiple regression

The scheme of our viscoelastic model is: vector e → {multiple regression Equation (13)} → viscoelastic coefficients v → {viscoelastic model Equation (4)} → moduli E ′ L and G ′ LT → {metamodel Equation (1)} → reconstructed resonance frequencies r .

It has a very interesting property: the same model of cane can be used for any kind of reeds (for instance bass clarinet or saxophone) or for any other boundary conditions. Only Equation (1) has to be changed (or at least, the coefficients a m,q have to be recomputed).

Within our particular setup, the viscoelastic model is however not required for reconstructing the observed resonance frequencies: PCA is theoretically the optimal linear scheme, in terms of least mean square error, for compressing a set of high dimensional vectors into a set of lower dimensional vectors and then reconstructing the original set by multiple regression. The shortened scheme is merely: vector f actor → {multiple regression} → reconstructed resonance frequencies r .

Let us examine this option. For this purpose we use the array f actor n,k computed in Section 3.2, holding our 4 principal components (factors). As before (Section D.3), we set f actor n,0 = 1. The array of reconstructed resonance frequencies r n,s,m can be computed by multiple regression using an array of matrices W s,m,k :

r[n, s] = W [s] • f actor[n]

(20)

As before (Section D.4), we have also the option to reduce the dimensionality from K to Ǩ, using the previously defined array of matrices M [s] and then use a unique matrix Ŵm, ǩ to operate the multiple regression:

r[n, s] = Ŵ • M [s] • f actor[n] (21)
We call these 2 options: hypotheses H10 and H11. For H11, we performed a small orthogonal rotation of the factors to concentrate the information about the hygrometric material properties in f actor n,4 , for a better fit. The results are summarized in Table 11.

Multiple regression is an accurate way to retrieve our measurements, comparable to H8 and H9. With only 48 coefficients, H11 is very efficient, even better for "transversal modes" as H10. The results of the regressive model are more difficult to interpret than those of the viscoelastic model. As e n,4 before, f actor n,4 serves uniquely to adjust f actor n,2 relatively to series B. As expected, f actor n,1 influences mainly the "transversal modes" and f actor n,2 the "flexural modes". Their respective coefficients in the matrix Ŵ reflect this antinomy. "Transversal" and "flexural modes" are concerned by f actor n,3 in a quite similar way, but the slope is not same.

Figure 1 :

 1 Figure 1: Holographic setup. L: main laser; AOM1, AOM2: acousto-optic modulators; M: mirror; BS: beam splitter; BE: beam expander; CCD: CCD camera; LS: loudspeaker exciting the clarinet reed through the bore of a clarinet mouthpiece at frequency f = ω/2π.

Figure 4 :

 4 Figure 4: (a) Flexural modes F 1, F 2, F 3, F 4 and F 5. Side view. (b) transversal modes T 1, T 2, T 3, T 4, T 5. Front view. (c) Generic modes: X1, X2, X3, X4, X5 and X6. View from above. The intersections of nodal lines with the sides of the reed are symbolized by the blue dots.

Figure 5 :

 5 Figure 5: The clarinet reed: coordinates system. (a) top view, (b) front view, (c) side view. x, y, z: Cartesian axes of the object. L, T, R: axes of the orthotropic material (L: longitudinal, T: tangential, R: radial).In the software Catia used for the simulations, the orthotropic model is Cartesian and not cylindrical. Therefore an exact equivalence between x, y, z, and L, T, R respectively can be assumed since we observed no important deviation between the direction of the grain and the axe of symmetry of the reed.

Figure 6 :

 6 Figure 6: Qualitative comparison with FEM computation, and typical variability of the experimental results. (a) First row: Quasi-static pattern (at 605 Hz, with strong excitation by the LS). 2nd to 5th row: Flexural modes F 1, F 2, F 3 and F 4. Leftmost column: numerical simulation of eigenmodes by FEM. Columns 2 to 7: Arbitrary selection representing the observed variability (The first two rows correspond to the same selection of reeds). Notice the marked asymmetries and the differences in the curvature of the interference fringes near the tip of the reed. (b) Rows 1 to 5: transversal modes T 1, T 2, T 3, T 4 and T 5 (probably). Columns: see (a). (c) Rows 1 to 6: Generic modes X1, X2, X3, X4, X5 and X6. The identification with X6 is sometimes unlikely. Columns: see (a).

Figure 8 :

 8 Figure 8: Schematic representation of the standard linear solid: two springs E 1 , E 2 and a dashpot E 3 .

  e 1 [n], e 2 [n], e 3 [n] and e 4 [n] are the 4 independent components characterizing the mechanical properties of the reed n. The choice of the notations is as follows: e k [n] is the kth element of the vector e[n]., which depends on n. These components are conditioned similarly to PCA as orthogonal factors: mean 0, standard deviation 1 and intercorrelation 0. The elimination of the moisture content can be achieved by reducing the components to a number of 3 for each series s = 1 (series A) and s = 2 (series B): for reasons explained in Appendix D.4, these components are denoted ě1 [s, n], ě2 [s, n] and ě3 [s, n] .

Figure 9 :

 9 Figure 9: Hypothesis H4: plot of the storage moduli E ′ L (ω) and G ′ LT (ω) [in MPa], according to Equations (4 to 7), computed for the mean value of all reeds. (a): E ′ L series A, (b): E ′ L series B, (c): G ′ LT series A and (d): G ′ LT series B. For G ′ LT , the moduli are multiplied by 10. The abscissa of the relaxation frequencies [in Hz] is denoted by dashed lines (series A) and dotted lines (series B). The corresponding numerical values are listed in Tables4 and 5. Only the portions of the curves between 2 and 18 kHz could be fitted adequately. The curves outside this range are purely hypothetical: we have no measurements.

d 22 =

 22 2d 12 ; d 13 = d 12 = -d 23 ; d 40 = d 50 = 6d 41 ; d 51 = 2d 41 = -d 53 ; (19) c 22 = c 24 = 1/2.

Table 1 :

 1

	1568, 2093, 2637, 3136,

Correlation (loadings)

between rotated factors from PCA and : i) variables (measured resonance frequencies) or ii) components from viscoelastic model (e[n], see Section 5.2), for comparison, sorted in reverse order of magnitude. In bold: greater correlation for each variable.

Table 2 :

 2 One-At-a-Time sensitivity study by FEM: Averaged ratio between relative change in frequency and relative change for each elastic coefficient (i.e. ±10%), sorted by decreasing order of magnitude, for the first 16 eigenmodes. F : flexural modes (F 1 to F 4), T : transversal modes (T 1 to T 4), X: generic modes (X1 to X6), L: lateral modes (L1 and L2; these modes were not observed in our study), All modes: averaged ratio over all modes. In bold: maximum absolute value for each coefficient of the orthotropic material: E L , E T and E

		F	T	X	L	All modes
	E L	0.4087	0.1053	0.1835 0.2093	0.2235
	G LT	0.0140	0.2681 0.1962 0.1067	0.1575
	E T	0.0076	0.0976 0.0741 0.0166	0.0562
	G LR	0.0438	0.0131	0.0257 0.0818	0.0341
	E R	0.0176	0.0120	0.0135 0.0662	0.0207
	G T R	0.0046	0.0092	0.0077 0.0215	0.0091
	ν T R	0.0015	0.0009	0.0015 0.0054	0.0019
	ν LT	0.0018 -0.0031 0.0004 0.0007	-0.0001
	ν LR	0.0009	0.0002	0.0005 0.0013	0.0007

R : Young moduli; ν LT , ν LR and ν T R : Poisson coefficients; G LT , G LR and G T R : shear moduli.

Table 3 :

 3 Optimal 9-coefficients model (hypothesis H4): Statistics of the frequency-dependent increase in storage modulus (for Series A or B and E L or G LT ), between the 1 st and the 4 th flexural mode: 1

	modulus	Series A Series B Series A Series B
	Mean	31%	16%	38%	36%
	Standard deviation	5%	6%	10%	9%
	Minimum	20%	1%	13%	15%
	Maximum	44%	32%	52%	49%

Table 4 :

 4 Summary of our results for hypotheses H4 and H5 about the viscoelastic behavior of the longitudinal Young modulus E L in cane. E 1 , E 2 and f r (relaxation frequency): parameters from Zener model. E ′ L at 4 kHz : storage modulus at 4 kHz [in MPa]. SD: standard deviation (the value preceding SD is the average). The model is valid only between 2 and 18 kHz.

	4 kHz

Table 5 :

 5 Summary of our results for hypotheses H4 and H5 about the viscoelastic behavior of the shear modulus in longitudinal / tangential plane G LT in cane. G ′ LT at 4 kHz: storage modulus at 4 kHz. Same structure as Table4. The model is valid only between 2 and 18 kHz.

	Model	Mode AF 1 AF 2 AF 3 Mode BF 1 BF 2 BF 3
		AF 2	0.73	BF 2	0.63
	observations	AF 3	0.34 0.75	BF 3	0.43	0.84
		AF 4	0.28 0.65 0.74	BF 4	0.00	0.43	0.45
		AF 2	0.93	BF 2	0.79
	H4	AF 3	0.69 0.90	BF 3	0.62	0.97
		AF 4	0.56 0.81 0.98	BF 4	0.36	0.82	0.91
		AF 2	0.90	BF 2	0.69
	H5	AF 3	0.57 0.87	BF 3	0.52	0.98
		AF 4	0.43 0.77 0.98	BF 4	0.26	0.84	0.92

Table 6 :

 6 Correlations of resonance frequencies between flexural modes within the same series, after observations and viscoelastic models H4 and H5. AF 1, AF 2, AF 3, BF 1, BF 2 and BF 3: correlations with modal frequencies F 1, F 2 and F 3 within series A or series B. Lines: corresponding flexural mode between which the correlations are computed.

Table 8 :

 8 Network of points for interpolating the thickness of the vamp. See explanations in the text.

	x [mm]	s 0 [mm] y=0 mm	s 1 [mm] y=4mm	s 2 [mm] y=6mm
	0	0.074	0.080	0.042
	5	0.343	0.293	0.197
	10	0.648	0.542	0.377
	15	1.047	0.847	0.571
	20	1.451	1.135	0.745
	25	1.926	1.527	1.078
	30	2.540	2.084	1.589
	35	3.351	2.817	2.256

Table 9 :

 9 Coefficients of the metamodel, Equation (1). The maximum negative and positive deviations of the model (compared to the values calculated by FEM) are given by δ-and δ+ [in cents]. The mode L1 (6th mode) has been deleted from theTable, since we didn't observe it.

	0 or x ≥ 67.5 x < 1.13196 -5.31 + 6.8x -x 2 x < 2.94661 (24.4 -x)x 4.08044 + √ 263 40 -11 900 x x < 67.5
	The thickness of the reed at point (x, y) is defined by:
	thickness(x, y) =	min[heel(y), vamp(x, y)] Abs(y) < contour(x) 0 otherwise
	C Coefficients of the metamodel

Table 10 :

 10 List of indices

Table 11 :

 11 2): 4 factors: 91.2% of the variance (53.6, 21.4, 10.8 and 5.4% for each factor). Synthesis of some hypotheses tested with our model. K: Number of components for each reed n; model: elastic, viscoelastic or multiple regression (with the elastic model, E L and G LT are independent from frequency; the elastic model is computed after the viscoelastic one by setting very small values for the coefficients affecting E 2 , E 3 , G 2 and G 3 , in order to avoid division by 0; multiple regression: see Equations (20 and 21)); #ActiveCoef : number of active coefficients in the matrices M or M and M (fitted through the fitting procedure); #OtherCoef : number of other coefficients estimated in model (a m,0 : adjusted so that the mean error for each mode m and both series is zero, W s,m,k : computed analytically, otherwise: fitted through the fitting procedure); RMSD [in cents] after Equation[START_REF] Obataya | Mechanical relaxation processes due to sugars in cane (Arundo donax L.)[END_REF]: this is a measure of goodness of fit (remember that the vibration patterns of the reeds were observed in steps by 25 cents).

	Hypothesis K	model	#ActiveCoef #OtherCoef RMSD
	H1 H2 H3 H4 H5 H6 H7 H8 H9 H10 H11	0 2 3 4 viscoelastic elastic elastic elastic 4 viscoelastic 4 viscoelastic 4 viscoelastic 4 viscoelastic 12 viscoelastic 4 regression 4 regression	2 4 7 9 18 44 60 60 --4	a m,0 → 11 a m,0 → 11 a m,0 → 11 a m,0 → 11 a m,0 → 11 a m,0 → 11 a m,0 → 11 a m,q → 66 a m,q → 66 W s,m,k → 110 Ŵm, ǩ → 44	76.2 54.8 43.9 30.4 29.8 29.1 28.6 23.2 19.8 21.8 24.7

The significantly lower correlations between resonance frequencies (compared to our data) shows that it was probably not the case in Pinard's study. This fact may also reflect an unprecise determination of the resonance frequencies.

The original notation from the cited paper is kept. This notation is only valid for this paragraph.

Do not confuse the morphological mode classes L1, L2, T 1, T 2, T 3 and T 4 with the axes L and T of the orthotropic material.

This notation allows to write the parameters of the model as a vector, as required by the computations, but unfortunately it hides the fact that the nature of E 3 (dashpot) is physically different from E 1 and E 2 (springs).

For this example, the component vector has 4 elements to be determined, plus a fifth element, which is a constant.

The parameter E 2 (respectively G 2 ) determines the influence of the Maxwell arm in the Zener model, since E 3 (respectively G 3 ) is constant in our simplified model; if E 2 = 0 the model is perfectly elastic and the viscous component disappears

Because H5 is probably less biased than H4 with its minimal structure

years. We also thank Bruno Gazengel, Morvan Ouisse and Emmanuel Foltête for fruitful discussions, and the referees for useful suggestions.

AF1/BF1 means AF1 versus BF1. The correlations are computed between AF 1[n] and BF 1[n], for n = 1 to N , where N = 55; missing observations are deleted.

This was not the case with PCA.

The fitting process was realized by repeatedly performing four kinds of procedures, in an arbitrary order :

These coefficients can be fitted through the fitting procedure (some constraints have however to be introduced, to avoid an important deviation from their theoretical values) or merely adjusted a posteriori, so that the total averaged deviation for each mode and both series is 0.

Acknowledgments

We express gratitude for the financial aid from association DEPHY of the "département de physique de l'ENS", and to Fadwa Joud for her help during the experiments. The Haute-Ecole ARC (Neuchâtel, Berne, Jura) provided the software for simulation and calculation, and offered helpful instructions. The Conservatoire neuchâtelois partially supported this study for nearly

APPENDICES

A Statistics and correlations

The table 7 gives the detailed statistics.

Linear correlations have been computed between all possible couples of variables in a usual way (all the variables from Table 7, except X5 and T 5 of series A, i.e. 13 variables for series A and 15 for series B, plus the nominal reed strength). Results of series A for the mode F 1 are denoted AF1, and similarly for the other results. Between the two series, the correlation is excellent for corresponding transversal modes (T 1 (i.e. AT1/BT1): 0.97, T 2: 0.97, T 3: 0.96 and T 4: 0.98), and generally good for corresponding generic modes (X1: 0.87, X2: 0.84, X3: 0.87 and X4: 0.55). For flexural modes, the correlation is good for F 1 and progressively lower for increasing mode order (F 1: 0.92, F 2: 0.66, F 3: 0.57 and F 4: 0.47).

Within the same series, on the contrary, there is a poor correlation between AF1 and all measurements of series A, and similarly for BF1 and series B. This is striking: the two best correlated variables are AF2 and AX1 (0.73 and 0.49, respectively) for AF1 and BF2 and BT1 (0.63 and 0.49, respectively) for BF1. Moreover, these correlations are quite low among all flexural modes: see Table 6. This fact is discussed in section 5.3.

The nominal reed strength correlates at 0.7 with AT1 and AX4. We expected a better correlation with F 1 (only 0.6). This is surprising, since the reeds were probably sorted by a quasi-static bending method by the manufacturer. This would mean that the storage modulus of E L at very low frequency is not well correlated with its value at the frequencies of the measured resonances. The influence of density has also to be considered. However, Obataya et al. [START_REF] Obataya | The effects of water soluble extractives on the acoustic properties of reed (Arundo donax L.)[END_REF] observed a good correlation between density and E L . This point has to be investigated (see also section D.4 in Appendix D).