# Characterisation of Anti-Arrhythmic Drug Effects on Cardiac
Electrophysiology using Physics-Informed Neural Networks

###### Abstract

The ability to accurately infer cardiac electrophysiological (EP) properties
is key to improving arrhythmia diagnosis and treatment. In this work, we
developed a physics-informed neural networks (PINNs) framework to predict how
different myocardial EP parameters are modulated by anti-arrhythmic drugs.
Using in vitro optical mapping images and the 3-channel Fenton-Karma model, we
estimated the changes in ionic channel conductance caused by these drugs.

Our framework successfully characterised the action of drugs HMR1556,
nifedipine and lidocaine – respectively, blockade of $I$K, $I$Ca, and $I$Na
currents – by estimating that they decreased the respective channel
conductance by $31.8\pm 2.7\%$ $(p=8.2\times 10^{-5})$, $80.9\pm 21.6\%$
$(p=0.02)$, and $8.6\pm 0.5\%$ $(p=0.03)$, leaving the conductance of other
channels unchanged. For carbenoxolone, whose main action is the blockade of
intercellular gap junctions, PINNs also successfully predicted no significant
changes $(p>0.09)$ in all ionic conductances.

Our results are an important step towards the deployment of PINNs for model
parameter estimation from experimental data, bringing this framework closer to
clinical or laboratory images analysis and for the personalisation of
mathematical models.

Index Terms—  Cardiac Electrophysiology, Physics-Informed Neural Networks
(PINNs), Parameter Identification, Optical Mapping, Mathematical Modelling

##  1 Introduction

Physics-Informed Neural Networks (PINNs) are a machine learning method that
combines knowledge of the equations of a system with data learning [1]. This
domain knowledge allows PINNs to learn from a small fraction of the data that
conventional neural networks require and ensures predictions consistent with
physics models. This makes PINNs a suitable candidate for biomedical
applications, where data are often sparse and noisy and interpretability is
highly desired.

Within cardiac EP, Sahli-Costabal et al. used PINNs to solve the isotropic
diffusion eikonal equation, using in silico data. They estimated arrival times
of the action potential (AP) and conduction velocity maps across the left
atrial surface [2]. Following this work, Grandits et al. solved the
anisotropic eikonal equation to learn fibre orientations and conductivity
tensors, using synthetic data and data from one patient [3]. Closely related
to our work, Herrero Martin et al. used PINNs with the monodomain equation on
sparse maps of transmembrane potential to estimate EP parameters, such as the
isotropic diffusion coefficient or surrogates of the AP duration [4]. The
study relied on the Aliev Panfilov model, a simple 2-variable model without
clear biological counterparts [5]. To provide clinically useful
characterisation of cardiac EP properties under the influence of drugs, PINNs
will need to be deployed in EP models that better capture ionic channel
complexity and demonstrate its ability to make inferences from experimental
data.

To this end, we coupled PINNs with the Fenton-Karma model, which is a
3-variable model of the cardiac AP with formulations of three transmembrane
currents: a fast inward $I_{\text{fi}}$, slow outward $I_{\text{so}}$, and
slow inward $I_{\text{si}}$ [6]. They are analogous to Na+, K+, and Ca2+
currents, respectively. We selected this model for its ability to reproduce
restitution properties and capture arrhythmic conditions, with the appropriate
level of complexity for PINNs. The Fenton-Karma equations are as follows:

$\partial_{t}u=\nabla\cdot(D\nabla
u)-J_{\text{fi}}(u;v)-J_{\text{so}}(u)-J_{\text{si}}(u;w),$ (1a)
$\partial_{t}v=\Theta(u_{c}-u)(1-v)/\tau_{v}^{-}(u)-\Theta(u-u_{c})v/\tau_{v}^{+},$
(1b)
$\partial_{t}w=\Theta(u_{c}-u)(1-w)/\tau_{w}^{-}-\Theta(u-u_{c})w/\tau_{w}^{+},$
(1c)

where the three scaled currents are given by

$J_{\text{fi}}(u,v)=-\frac{v}{\tau_{d}}\Theta(u-u_{c})(1-u)(u-u_{c}),$ (2a)
$J_{\text{so}}(u)=\frac{u}{\tau_{o}}\Theta(u_{c}-u)+\frac{1}{\tau_{r}}\Theta(u-u_{c}),$
(2b)
$J_{\text{si}}(u;w)=-\frac{w}{2\tau_{\text{si}}}(1+\tanh[k(u-u_{c}^{\text{si}})]).$
(2c)

Here, $u$ is the dimensionless electrical potential across the cell membrane,
scaled to take values from $0$ to $1$. $v$ and $w$ are the latent gate
variables for $I_{\text{fi}}$ and $I_{\text{si}}$, respectively. $\Theta(x)$
is the Heaviside step function defined by $\Theta(x)=1$ for $x\geq 0$ and
$\Theta(x)=0$ for $x<0$. Membrane potential diffuses to the neighbouring
cells, as modelled by the diffusion term for $u$, $\nabla(D\nabla u)$. We
assume homogeneous and isotropic conduction, in which case $D$ is a constant
scalar. The 8 $\tau$ parameters are various time constants. Among them,
$\tau_{d}$, $\tau_{r}$, $\tau_{o}$, and $\tau_{si}$ are inversely proportional
to ionic channel conductances, as shown in Eq. (2a)-(2c). We model the effects
of anti-arrhythmic drugs as time-independent reductions in channel
conductances, that is, as increases in the $\tau$ variables.

Aims. We aim to demonstrate how PINNs can characterise the effects of anti-
arrhythmic drugs on ionic channels using in vitro optical mappings images from
cardiomyocyte preparations. We will use PINNs to solve the Fenton-Karma
equations and simultaneously estimate model parameters related to ionic
channel conductances. For this we will use:

  1. 1.

in silico data, which allow us to assess the accuracy and tune the model.

  2. 2.

in vitro optical mapping data from cardiomyocyte preparations.

##  2 Methods

All code used in this study is available at github.com/annien094/EP-PINNs-for-
drugs.

###  2.1 In Silico Data Generation

We generated in silico membrane potentials on a 1D cable of $10$mm by solving
the Fenton-Karma model with a central finite difference and an explicit
4-stage Runge-Kutta method, using a time step of $5$ and a spatial step of
$100$. The values of model parameters were set as in [7] (see FK-CAF model).
We used APs across time from two spatial points as inputs to the PINNs, to
mimic the available experimental data. We assessed the accuracy of parameter
estimation using the relative error (RE).

###  2.2 Optical Mapping Images

We tested the performance of our PINNs setup on in vitro optical mapping
images from neonatal rat ventricular myocyte preparations. These preparations
were stained with a voltage-sensitive dye (di-8-ANNEPS) before imaging at a
high spatio-temporal resolution of $400$$\times$$85$ pixels at $525.39$
frames/sec to obtain uncalibrated measurements of electrical potentials. All
experimental details are described in [8].

We investigated the effects of four drugs: HMR1556, nifedipine, lidocaine, and
carbenoxolone (CBX). Their effects are summarised in Table 1. We assessed the
effects of drugs on PINNs’ estimates of model parameters. We identified drug-
induced significant differences using a one-tailed two-sample $t$-test at a
significance level of $0.05$.

| HMR1556 | Nifedipine | Lidocaine | CBX
---|---|---|---|---
Site of blockade |  Potassium channel ($I$K) |  Calcium channel ($I$Ca) |  Sodium channel ($I$Na) |  Intercellular gap junctions
Effects on model parameters |  Increases $\tau_{r}$ |  Increases $\tau_{si}$ |  Increases $\tau_{d}$ |  No effect on all $\tau$
Effects on AP |  Prolongs duration |  Shortens duration |  Prolongs upstroke |  Decreases conduction velocity
Table 1: Effects of anti-arrhythmic drugs [9, 10, 11, 12].

For each drug, we used four time series of optical mapping images. Two of
these series were acquired under the action of the drug at half-maximal
inhibitory concentration (IC50), at two different preparation locations. The
other pair was acquired in the absence of drugs (baseline conditions) also at
two different locations. For each time series, we improved the signal-to-noise
ratio of the data by:

  * •

spatially averaging the signal over a circular area of radius $50$ (indicated
in Fig. 1);

  * •

applying a temporal moving average filter to the signal except at the AP
upstrokes;

  * •

correcting temporal baseline drifts by subtracting each AP by the average of
$10$ before each upstroke;

  * •

overlaying and averaging all APs across time;

  * •

and scaling each averaged AP to the [0, 1] interval.

Fig. 1: Data processing pipeline for optical mapping images.

For simplicity, we modelled the cardiomyocyte preparation as a $10$mm cable
with the measured APs placed at positions $x_{1}=2$mm and $x_{2}=8$mm. A
diffusion coefficient of $D=0.1$ was used. These were then used as inputs for
PINNs. This pipeline is illustrated in Fig. 1.

###  2.3 PINNs Setup

We used a fully-connected neural network with 4 layers of 32 neurons each, a
$tanh$ activation, a Glorot Uniform initialisation, and Adam optimisation.
There were two input variables, space and time position $x$ and $t$, and three
output variables: $u$, $v$, and $w$. We trained the network for $100,000$
epochs with learning rate $5\times 10^{-5}$. These hyperparameters were
selected empirically based on tests using in silico data. The loss function
consists of 8 terms:

$\displaystyle L$ $\displaystyle=L_{fu}+L_{fv}+L_{fw}+L_{\text{BC}}+$
$\displaystyle
L_{\text{IC},u}+L_{\text{IC},v}+L_{\text{IC},w}+L_{\text{data}}.$ (3)

$L_{fu}+L_{fv}+L_{fw}$ accounts for the agreement with the Fenton-Karma
equations Eq. (1a)-(1c). $L_{\text{BC}}$ and
$L_{\text{IC},u}+L_{\text{IC},v}+L_{\text{IC},w}$ accounts for the no-flux
Neumann boundary condition for membrane voltage $\frac{\partial
u}{\partial\vec{n}}=0$ and initial conditions for variables $u$, $v$, and $w$.
$L_{\text{data}}$ accounts for the agreement with experimental measurements of
$u$. We empirically determined that equal weights for all terms leads to the
best network performance.

Fig. 2: The PINNs architecture. Inputs $x$ and $t$ are the spatial and time
points. The network optimises its weights by minimising the loss function $L$,
to provide solutions for $u$, $v$, and $w$. It simultaneously estimates model
parameter(s) $\tau$.

For training, we used $60000$ points within the domain for the agreement with
equations, $1000$ boundary points, and $99$ initial points for boundary and
initial conditions. We split all available data into train and test sets at
$9$:$1$ and calculate the root mean squared error (RMSE) on the test points.

###  2.4 Model Parameter Estimation using PINNs

We used PINNs to estimate model parameters $\tau_{d}$, $\tau_{r}$, $\tau_{o}$,
and $\tau_{si}$ in two modes: a) individually, with the other parameters fixed
at literature values (see FK-CAF in [7]), and b) estimating all four
parameters simultaneously. We initialised the model parameter(s) by randomly
drawing from a uniform distribution around the ground truth(s), where the
distribution widths were based on uncertainties from preliminary results.

To estimate the uncertainty in parameter estimates, the network was randomly
initialised to a different configuration five times and the standard deviation
of the model parameter estimates over these runs was computed.

##  3 Results

###  3.1 Parameter estimation with in silico data

Fig. 3: PINNs’ predictions and in silico ground truth data for transmembrane
potential $u$ and gate variables $v$ and $w$ at $2$ on a 1D $10$ cable. Only
data for $u$ were given here and at $8$, and none for $v$ and $w$. Bottom
right: 10 independent runs of parameter estimates of $\tau_{d}$ across epochs.

Using PINNs, we were able to accurately estimate the Fenton-Karma model
parameters, and accurately reproduce the APs and time course of latent
variables, as shown in Fig. 3. PINNs predictions for variables $u$, $v$, and
$w$ are in good agreement with the synthetic ground truth data, even when only
data for $u$ were given at two locations. Estimates of the variables were not
affected by the number of model parameters predicted by PINNs: RMSE was
$(4.6\pm 2.0)\times 10^{-3}$ when estimating each parameter individually, and
$(4.9\pm 2.2)\times 10^{-3}$ when estimating all four simultaneously.

As shown in Fig. 4, when estimating each parameter separately, REs were within
$15\%$. REs were generally larger when estimating four parameters
simultaneously, as expected, but even in this regime did not exceed $30\%$.
$\tau_{si}$ (inverse of $I$Ca conductance) was the least accurate in both
modes, and $\tau_{d}$ (inverse of $I$Na conductance) the most accurate.

Fig. 4: Relative errors of parameter estimates using in silico data.

###  3.2 Drug Effect characterisation using PINNs

With experimental in vitro data as inputs, PINNs were also able to accurately
reproduce the experimental APs, with an average RMSE for $u$ within $7\times
10^{-2}$ (Fig. 5).

Our method remarkably identified the effects of drugs on each ionic channel by
measuring significant increases in the $\tau$ parameters that correspond to
the ionic channels each drug is known to modulate. The relative changes in
$\tau$ and the associated $p$-values are reported in Table 2. We estimated
each parameter separately for in vitro data for stability.

Fig. 5: PINNs’ predictions of the action potentials and in vitro optical
mapping data showed good agreement even when some model parameters were not
given to the network. The RMSE for $u$ in baseline conditions was $(2.6\pm
2.5)\times 10^{-2}$.

| HMR1556 | Nifedipine
---|---|---
$\Delta\tau_{d}$ | $-3.4\pm 0.3\%$ | $-5.5\pm 0.5\%$
$p$-value | 0.50 | 0.28
$\Delta\tau_{r}$ | $31.8\pm 2.7\%$ | $-0.69\pm 0.04\%$
$p$-value | $8.2\times 10^{-5}$ | 0.85
$\Delta\tau_{o}$ | $37.1\pm 14.5\%$ | $9.5\pm 1.4\%$
$p$-value | 0.065 | 0.24
$\Delta\tau_{si}$ | $-13.6\pm 2.6\%$ | $80.9\pm 21.6\%$
$p$-value | 0.19 | 0.024
Expected |  Increases $\tau_{r}$ |  Increases $\tau_{si}$

| Lidocaine | CBX
---|---|---
$\Delta\tau_{d}$ | $8.6\pm 0.5\%$ | $-8.3\pm 4.6\%$
$p$-value | 0.03 | 0.73
$\Delta\tau_{r}$ | $1.2\pm 0.2\%$ | $9.1\pm 0.9\%$
$p$-value | 0.84 | 0.09
$\Delta\tau_{o}$ | $17.7\pm 10.6\%$ | $-8.1\pm 6.2\%$
$p$-value | 0.58 | 0.81
$\Delta\tau_{si}$ | $-9.8\pm 2.0\%$ | $-2.2\pm 0.2\%$
$p$-value | 0.22 | 0.56
Expected |  Increases $\tau_{d}$ |  No change in $\tau$

Table 2: Relative changes in the $\tau$ parameters (inverse of channel
conductance) between baseline and drug and associated $p$-values. Significant
increases in $\tau$ are shown in green.

##  4 Discussion

In this work, we developed a PINNs framework capable of accurate
characterisation of anti-arrhythmic drug actions on EP parameters, using
experimental in vitro images. We coupled PINNs with the Fenton-Karma model,
which distinguishes among currents of the three main ion species in the heart:
Na+, K+, and Ca2+. Starting with in silico data, PINNs accurately estimated
the model parameters corresponding to different ionic channel conductances.
Then, with in vitro optical mapping images, PINNs identified significant
increases in model parameters corresponding to the expected drug-induced
reduction of channel conductances. In both cases, PINNs gave excellent
predictions for AP propagation, which demonstrated its capability to work well
with sparse measurements, experimental noise and artefacts.

So far, studies have mostly focussed on demonstrating PINNs performance on in
silico data [2, 13, 14]. We showed that PINNs are also effective for model
calibration with in vitro data as inputs. Previous work had used PINNs to
characterise drug effects with a simple model that does not consider ionic
channels [4]. Our work shows that it is possible to use more sophisticated
cardiac EP models to gain more detailed insights into drug mechanisms, which
has important applications in the safety and efficacy testing of new drugs
[15]. Compared to traditional numerical methods, our approach has the
advantages of being mesh-free, using automatic differentiation, and can be
compactly implemented to simultaneously solve forward and inverse problems
[14].

Limitations in our current framework include experiment-specific results,
where evaluation on other sets of data would be valuable for robustness. In
addition, the parameter estimates and convergence properties are sensitive to
the initialisation. Future work may focus on exploring alternative
optimisation schemes and loss function compositions. We also plan to extend
this framework to the analysis of electrograms to bring PINNs’ advantages
closer to the clinical setting.

##  5 Compliance with ethical standards

This is a computational study for which no ethical approval was required. The
data used were acquired in [8], with the appropriate ethical approvals
detailed therein: all procedures were conducted according to the standards set
by the EU Directive 2010/63/EU and were approved by the Imperial College
London Ethical Committee.

##  6 Acknowledgments

This work was supported by the British Heart Foundation (RE/18/4/34215,
RG/16/3/32175, PG/16/17/32069 and Centre of Research Excellence), the National
Institute for Health Research (UK) Biomedical Research Centre and the
Rosetrees Trust through the interdisciplinary award ”Atrial Fibrillation: A
Major Clinical Challenge”.

## References

  * [1] M. Raissi, P. Perdikaris, and G. E. Karniadakis,  “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,”  Journal of Computational Physics, vol. 378, pp. 686–707, 2 2019.
  * [2] Francisco Sahli Costabal, Yibo Yang, Paris Perdikaris, Daniel E Hurtado, and Ellen Kuhl,  “Physics-informed neural networks for cardiac activation mapping,”  Frontiers in Physics, vol. 8, pp. 42, 2020.
  * [3] Thomas Grandits, Simone Pezzuto, Francisco Sahli Costabal, Paris Perdikaris, Thomas Pock, Gernot Plank, and Rolf Krause,  “Learning atrial fiber orientations and conductivity tensors from intracardiac maps using physics-informed neural networks,”  in International Conference on Functional Imaging and Modeling of the Heart. Springer, 2021, pp. 650–658.
  * [4] Clara Herrero Martin, Alon Oved, Rasheda A Chowdhury, Elisabeth Ullmann, Nicholas S Peters, Anil A Bharath, and Marta Varela,  “Ep-pinns: Cardiac electrophysiology characterisation using physics-informed neural networks,”  Frontiers in Cardiovascular Medicine, vol. 8, pp. 768419, 2022.
  * [5] Rubin R Aliev and Alexander V Panfilov,  “A simple two-variable model of cardiac excitation,”  Chaos, Solitons & Fractals, vol. 7, no. 3, pp. 293–301, 1996.
  * [6] Flavio Fenton and Alain Karma,  “Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation,”  Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 8, no. 1, pp. 20–47, 1998.
  * [7] Amy M Goodman, Robert A Oliver, Craig S Henriquez, and Patrick D Wolf,  “A membrane model of electrically remodelled atrial myocardium derived from in vivo measurements,”  EP Europace, vol. 7, no. s2, pp. S135–S145, 2005.
  * [8] Rasheda A Chowdhury, Konstantinos N. Tzortzis, Emmanuel Dupont, Shaun Selvadurai, Filippo Perbellini, Chris D. Cantwell, Fu Siong Ng, Andre R. Simon, Cesare M. Terracciano, and Nicholas S. Peters,  “Concurrent micro-to macro-cardiac electrophysiology in myocyte cultures and human heart slices,”  Scientific Reports, vol. 8, no. 1, pp. 1–13, 12 2018.
  * [9] George P Thomas, Uwe Gerlach, and Charles Antzelevitch,  “Hmr 1556, a potent and selective blocker of slowly activating delayed rectifier potassium current,”  Journal of cardiovascular pharmacology, vol. 41, no. 1, pp. 140–147, 2003.
  * [10] Théophile Godfraind,  “Discovery and development of calcium channel blockers,”  Frontiers in pharmacology, vol. 8, pp. 286, 2017.
  * [11] BRUCE P Bean, CHARLES J Cohen, and RICHARD W Tsien,  “Lidocaine block of cardiac sodium channels,”  The Journal of general physiology, vol. 81, no. 5, pp. 613, 1983.
  * [12] Pipin Kojodjojo, Prapa Kanagaratnam, Oliver R Segal, Wajid Hussain, and Nicholas S Peters,  “The effects of carbenoxolone on human myocardial conduction: a tool to investigate the role of gap junctional uncoupling in human arrhythmogenesis,”  Journal of the American College of Cardiology, vol. 48, no. 6, pp. 1242–1249, 2006.
  * [13] Alireza Yazdani, Lu Lu, Maziar Raissi, and George Em Karniadakis,  “Systems biology informed deep learning for inferring parameters and hidden dynamics,”  PLoS Computational Biology, vol. 16, no. 11, pp. 1–19, 2020.
  * [14] Lu Lu, Xuhui Meng, Zhiping Mao, and G. E. Karniadakis,  “DeepXDE: A deep learning library for solving differential equations,”  SIAM Review, vol. 63, no. 1, pp. 208–228, 2 2021.
  * [15] Nejib Zemzemi, Miguel O Bernabeu, Javier Saiz, Jonathan Cooper, Pras Pathmanathan, Gary R Mirams, Joe Pitt-Francis, and Blanca Rodriguez,  “Computational assessment of drug-induced effects on the electrocardiogram: from ion channel to body surface potentials,”  British journal of pharmacology, vol. 168, no. 3, pp. 718–733, 2013.


# A method to estimate the rectangular orthotropic plate elastic constants
using least-squares and Chladni patterns

Michele Ducceschi  Department of Industrial Engineering, University of
Bologna, Italy  Sebastian Duran  Department of Industrial Engineering,
University of Bologna, Italy  Henna Tahvanainen  Department of Industrial
Engineering, University of Bologna, Italy  Ludovico Ausiello  SENE, University
of Portsmouth, UK

###### Abstract

A method to retrieve the elastic constants of rectangular wooden plates is
presented, relying on the measurement of a set of eigenfrequencies and the
identification of the corresponding mode shapes, and belonging to the more
general category of non-destructive inverse parameter estimation methods.
Compared to previous works, the current method is effective with any choice of
boundary conditions. Furthermore, the error function is linear in the elastic
constants, which are computed via a matrix inversion. This framework lends
itself naturally to a physical interpretation of the results in terms of
linear combinations of eigenmodes, yielding new sets of modes and associated
combined mode shapes in which the elastic constants are completely uncoupled.
A number of numerical benchmark tests and experimental cases are treated in
detail, highlighting the reliability of the proposed methodology in cases of
interest in acoustics and musical acoustics.

##  1 Introduction

.

Inverse parameter estimation refers to a wealth of techniques used to identify
model parameters in materials such as wood [1, 2], composites [3, 4], ceramics
[5, 6], metals [7, 8] and others. Commonly, one looks for parameters
minimising the discrepancy between the experimental data and a suitable
reference model, using various techniques differing in application and
methodology [9]. A typical distinction in this sense considers destructive and
non-destructive techniques. The measurement of the elasticity constants
through traction, compression, bending and twisting is a typical example of
the former: for elastic specimens, some form of permanent damage usually
results when the applied forces operate beyond the linear regime [8, 10].
Dynamic or vibrational testing, on the other hand, is an example of the
latter. In this case, the operating regime is most often linear, and the
measurement setup does not usually produce damage to the specimens [1, 2, 3,
4].

Inverse parameter estimation is useful to measure the parameter variability
across wood samples, the subject of this work. Wood is often modelled as an
orthotropic material, depending on nine independent elastic constants: three
Young’s moduli, three shear moduli, and three Poisson’s ratios [11]. Inverse
parameter techniques have been successfully devised for this case, for example
employing 3D laser vibrometry and finite element model updating [3, 12].
Differences across cuts from the same tree have been reported [13], [14,
ch.3], and [15, ch.4]. Furthermore, wood is susceptible to changes in
environmental conditions such as temperature and moisture [16], as well as to
chemical treatment [17, 18]. These changes entail macroscopic modifications in
a sample’s density and size, as well as in its response to excitation [19].
Whilst the former are easily recognised and measurable, the latter is
established by changes in the sample’s elastic properties, which are harder to
quantify.

In some cases, a reduction of the number of elastic constants is possible
through further approximations, such as when the sample is thin, a case most
often encountered in musical instruments where typical thicknesses are found
in the range of a few millimetres [20]. Well-established techniques for
obtaining the elastic constants of thin rectangular plates exist, such as the
ones presented by McIntyre and Woodhouse [1] and by Caldersmith and Freeman
[2]. In these cases, the elastic constants are retrieved by measuring the
frequencies of three specific modes of a plate with free edges, and by
comparing them with approximate analytic solutions obtained via the Rayleigh
method. These techniques have a simple and inexpensive setup and have become
commonplace among scientists and instrument builders alike [21, 22, 23].

In this work, an extension of such techniques is offered, not relying on a
specific set of boundary conditions or a fixed set of modes, and leveraging
numerical simulation. It will be shown that, in the thin-plate approximation,
the modal frequencies depend linearly on the elastic constants, where the
linear coefficients can be extracted from a batch of numerical training plates
sharing the same boundary conditions and aspect ratio as the experimental
plate.

Once such modal coefficients are known, the elastic constants can be estimated
immediately via a matrix inversion in a least-square sense. Compared to well-
established techniques in the literature, the current method offers the
possibility to perform multiple estimates of the elastic constants, using
various combinations of boundary conditions and modes. Furthermore, the
mathematical problem is linear, thus avoiding the need for nonlinear solvers
or iterative methods for which the existence and uniqueness of the solution
may not be guaranteed [24]. One only needs to associate the experimental
frequency to a matching set of modal coefficients depending on the nodal lines
on the plate’s surface, which can be achieved quickly using the known method
by Chladni [25, 26]. Numerical benchmark tests and experimental results show
the applicability of the method to the measurement of the elastic constants of
spruce boards. The paper is structured as follows: Section 2 presents the
methodology, with a specific focus on the linear dependence of the non-
dimensional squared plate frequencies on the elastic constants. Section 3
presents a few numerical benchmark tests, highlighting the feasibility of the
proposed methodology under two different choices of boundary conditions.
Section 4 describes the experimental setup to obtain the frequencies of a
cantilever Finnish spruce plate and a clamped red spruce plate. Finally,
Section 5 presents the estimates of the elastic constants for the two
experimental plates and a discussion on the validity of the method.

##  2 Methodology

The flexural vibration of a thin, rectangular wood panel can be described by
the orthotropic Kirchhoff plate equation [27]. This is:

$\rho
L_{z}\,{\partial_{t}^{2}u}=-\left(D_{x}({\partial_{x}^{4}}+\nu_{yx}{\partial_{x}^{2}\partial_{y}^{2}})+D_{y}({\partial_{y}^{4}}+\nu_{xy}{\partial_{x}^{2}\partial_{y}^{2}})+D_{s}\,{\partial_{x}^{2}\partial_{y}^{2}}\right)u.$
(1)

Here, $u=u({\bf
x},t):\mathcal{V}\times\mathbb{R}^{+}_{0}\rightarrow\mathbb{R}$ is the
flexural displacement of the plate, a function of the spatial coordinates
${\bf x}\in\mathcal{V}=\\{(x,y)\,|\,0\leq x\leq L_{x},\,0\leq y\leq L_{y}\\}$,
as well as time $t\geq 0$. In the above, $\partial_{r}^{l}$ denotes the
$l^{th}$ partial derivative with respect to $r$. Boundary conditions of
classic type are imposed at the plate’s edges, such as free, clamped or
simply-supported. Let $\mathcal{B}$ denote the set of boundary conditions
imposed along the plate’s edges, such that, for example,
$\mathcal{B}=\\{\text{C-S-F-C}\\}$ denotes a plate with a clamped, a simply-
supported, a free and another clamped edge. For an edge perpendicular to the
$y$ direction, these are defined as follows [27]:

C: $\displaystyle u=\partial_{x}u=0,$ (2a) S: $\displaystyle
u=(\partial^{2}_{x}+\nu_{yx}\partial^{2}_{y})u=0,$ (2b) F:
$\displaystyle(D_{x}\partial^{3}_{x}+(D_{x}\nu_{yx}+D_{s})\partial_{x}\partial^{2}_{y})u=(\partial^{2}_{x}+\nu_{yx}\partial^{2}_{y})u=0,$
(2c)

with analogous definitions holding for an edge perpendicular to the $x$
direction (one only needs to swap $x$ and $y$ in the above). A further
condition must be imposed at a corner of two free edges, namely $u_{xy}=0$
(this is known as the Kirchhoff corner condition).

In the model above, constants appear as: $\rho$, the volume density of the
plate (in kg $\cdot$ m-3); the plate dimensions $L_{x},L_{y},L_{z}$ (in m, and
where the thickness $L_{z}$ is assumed constant across $\mathcal{V}$); the
three rigidity constants $D_{x}$, $D_{y}$, $D_{s}$ (in N $\cdot$ m). These are
defined by:

$D_{x}=\frac{E_{x}L_{z}^{3}}{12(1-\nu_{xy}\nu_{yx})},\qquad
D_{y}=\frac{E_{y}L_{z}^{3}}{12(1-\nu_{xy}\nu_{yx})},\qquad
D_{s}=\frac{G_{xy}L_{z}^{3}}{3},$ (3)

where $E_{x}$, $E_{y}$ are Young’s moduli (in Pa), $\nu_{xy}$, $\nu_{yx}$ are
dimensionless Poisson’s ratios, $G_{xy}$ is a shear modulus (in Pa). The
rigidity constants are defined in terms of five elastic constants ($E_{x}$,
$E_{y}$, $G_{xy}$, $\nu_{xy}$, $\nu_{yx}$). However, due to the symmetry of
the compliance matrix in the elasticity equations, only four of them are
independent [28]. Usually, one of the two Poisson’s ratios is set according
to:

$E_{x}\nu_{yx}=E_{y}\nu_{xy}.$ (4)

The wood’s fibre orientation is responsible for the orthotropic character of
lumber and results from the particular cutting technique employed. In quarter
sawing, logs are quartered lengthwise with the annual rings placed almost
perpendicular to the board faces [29]. This technique, whilst wasteful,
produces boards with an increased resistance against environmental factors
such as moisture and is most prized in musical instrument making.
Conveniently, the axes’ orientation may be referred back to the fibres. In the
following, $x$ will denote the direction along the grain (longitudinal); $y$
will denote the direction across the grain (radial); and $z$ will denote the
tangential direction; see also Figure 1. The stiffness in the longitudinal
direction is usually an order of magnitude higher than in the other two
directions, resulting in increased wave propagation velocity at reference
wavelengths.

Figure 1: Fiber direction convention in this article. $x$: longitudinal; $y$:
radial; $z$: tangential.

###  2.1 Non-dimensional model

Equation (1) depends on six constants: $\rho$, $L_{z}$, $E_{x}$, $E_{y}$,
$G_{xy}$, $\nu_{xy}$, as well as the side lengths $L_{x},L_{y}$. Inverting a
system with such a large number of parameters may be impractical and lead to
non-uniqueness of the resulting solution, an issue often emerging in inverse
parameter estimation methods [30]. In some cases, typically in machine
learning and, particularly, deep learning frameworks, interpretability is
unnecessary so long as the expected system’s output is produced [31]. In other
cases, such as here, a physical interpretation of the model may be obtained by
recasting the problem in a form needing a minimal amount of independent
parameters, as one may do immediately through non-dimensionalisation. To that
end, divide both sides of (1) by $\rho L_{z}$, and collect $D_{s}$ on the
right-hand side. This gives:

${\partial_{t}^{2}u}=-\kappa^{2}\left(({\partial_{x}^{4}}+\nu_{yx}\partial_{x}^{2}\partial_{y}^{2})\,p+({\partial_{y}^{4}}+\nu_{yx}\partial_{x}^{2}\partial_{y}^{2})\,q+\,{\partial_{x}^{2}\partial_{y}^{2}}\right)u,$
(5)

where $\kappa^{2}:=\rho^{-1}L_{z}^{-1}D_{s}$, $p:=D^{-1}_{s}D_{x}$,
$q:=D^{-1}_{s}D_{y}$. Non-dimensionalisation of the system proceeds as
follows. Define:

$u={A}^{\tfrac{1}{2}}\,\bar{u},\,\,x={A}^{\tfrac{1}{2}}\,\bar{x},\,\,y={A}^{\tfrac{1}{2}}\,\bar{y},\,\,t=\kappa^{-1}{A\bar{t}},$
(6)

where the overbar notation denotes a non-dimensional variable, and where
$A:=L_{x}L_{y}$ is the plate surface area. Thus, (5) becomes:

${\partial_{\bar{t}}^{2}\bar{u}}=-\left(({\partial_{\bar{x}}^{4}}+\nu_{yx}\partial_{\bar{x}}^{2}\partial_{\bar{y}}^{2})\,p+({\partial_{\bar{y}}^{4}}+\nu_{yx}\partial_{\bar{x}}^{2}\partial_{\bar{y}}^{2})\,q+\,{\partial_{\bar{x}}^{2}\partial_{\bar{y}}^{2}}\right)\bar{u},$
(7)

defined over the non-dimensional domain
$\mathcal{\bar{V}}=\\{(\bar{x},\bar{y})\,|\,0\leq\bar{x}\leq\sigma^{\tfrac{1}{2}},\,0\leq\bar{y}\leq\sigma^{-\tfrac{1}{2}}\\}$,
with $\sigma:=L_{y}^{-1}L_{x}$ denoting the aspect ratio of the plate.
Equation (7) depends only on the two ratios $p,q$, as well as $\sigma$.

Transforming (7) in the frequency domain (that is, using
$\partial_{\bar{t}}\rightarrow j\bar{\omega}$, $\partial_{\bar{x}}\rightarrow
j\bar{\gamma}_{\bar{x}}$, $\partial_{\bar{y}}\rightarrow
j\bar{\gamma}_{\bar{y}}$ for temporal frequency $\bar{\omega}$ and wavenumbers
$\bar{\gamma}_{\bar{x}}$, $\bar{\gamma}_{\bar{y}}$), the non-dimensional
dispersion relation is recovered:

$\bar{\omega}^{2}=(\bar{\gamma}_{\bar{x}}^{4}+\nu_{yx}\bar{\gamma}_{\bar{x}}^{2}\bar{\gamma}_{\bar{y}}^{2})\,p+(\bar{\gamma}_{\bar{y}}^{4}+\nu_{xy}\bar{\gamma}_{\bar{x}}^{2}\bar{\gamma}_{\bar{y}}^{2})\,q+\bar{\gamma}_{\bar{x}}^{2}\bar{\gamma}_{\bar{y}}^{2},$
(8)

from which the modal frequencies $\bar{\omega}_{m,n}$, are extracted as a
quantised version of (8). Here, $m,n\in\mathbb{N}$ denote a pair of modal
indices related to the nodal lines in the $x$ and $y$ directions. For fixed
$p$, $q$, the quantisation depends exclusively on the boundary conditions
$\mathcal{B}$ and the aspect ratio $\sigma$. On the other hand, when one
allows $p$, $q$ to vary keeping $\sigma$ and $\mathcal{B}$ fixed, a linear
dependence is assumed:

$\bar{\omega}^{2}_{m,n}=a_{m,n}\,p+b_{m,n}\,q+c_{m,n},$ (9)

for modal number-dependent coefficients $(a_{m,n},b_{m,n},c_{m,n})\geq 0$.
This is a key feature of the current formulation, allowing to estimate the
rigidity constants using a simple matrix inversion, as will be shown in the
following sections. Before proceeding, notice that the dimensional modal
frequencies $\omega_{m,n}$ are obtained via multiplication of the non-
dimensional frequencies by the time scaling constant, as per (6):

$\omega^{2}_{m,n}:=\kappa^{2}A^{-2}\bar{\omega}^{2}_{m}=D_{s}\,\rho^{-1}L_{z}^{-1}A^{-2}\bar{\omega}^{2}_{m,n}.$
(10)

Substituting the linear relationship (9) in (10), one obtains:

$\omega^{2}_{m,n}=\lambda\left(a_{m,n}D_{x}+b_{m,n}D_{y}+c_{m,n}D_{s}\right),$
(11)

with $\lambda:=\rho^{-1}L_{z}^{-1}A^{-2}$.

###  2.2 Problem Formulation

In the following, assume one wishes to estimate the rigidity constants
$D_{x}$, $D_{y}$ $D_{s}$ for the experimental plate under study. Furthermore,
the plate parameters $\rho$, $L_{x},L_{y},L_{z}$ are assumed to be known,
along with the boundary conditions $\mathcal{B}$.

The method relies on minimising a suitably defined error between the
dimensional numerical eigenfrequencies $\omega_{m,n}$ in (11) and a set of
experimentally measured frequencies $\hat{\omega}_{m,n}$, where the
minimisation is performed with respect to the unknown parameters
$d_{x},d_{y},d_{s}$ (these are approximations of the “true” rigidity constants
$D_{x}$, $D_{y}$, $D_{s}$). Thus, define:

${\boldsymbol{\psi}}:=\lambda^{-1}\hat{\boldsymbol{\omega}}^{2},\quad{\boldsymbol{\phi}}\left(d_{x},d_{y},d_{s}\right):=\,\begin{bmatrix}{\bf
a},{\bf b},{\bf c}\end{bmatrix}\begin{bmatrix}d_{x}\\\ d_{y}\\\
d_{s}\end{bmatrix}:={\bf X}\,{\bf d},$ (12)

where ${\bf X:=\begin{bmatrix}{\bf a},{\bf b},{\bf c}\end{bmatrix}}$, ${\bf
d:=\begin{bmatrix}d_{x},d_{y},d_{s}\end{bmatrix}^{\intercal}}$. Here, $\bf a$
is a column vector obtained by stacking $a_{m,n}$ for a set of modes. As an
example, assume to work with modes $(0,1),(1,1),(1,2)$. Then, ${\bf
a}=[a_{0,1},a_{1,1},a_{1,2}]^{\intercal}$. Analogous definitions hold for $\bf
b$, $\bf c$. Thus, $\boldsymbol{\phi}$ is a vector expression of (11) using
the approximate rigidity constants $\bf d$.

The error between the measured and the numerical frequencies is here defined
as:

$\varepsilon\left(d_{x},d_{y},d_{s}\right):=\frac{1}{2}\|{\boldsymbol{\psi}}-{\boldsymbol{\phi}}\|^{2}_{2},$
(13)

from which the rigidity constant vector $\bf d$ is obtained through a least-
square minimisation [32]:

${\bf d}=\arg\min\varepsilon=\left({\bf X}^{\intercal}{\bf
X}\right)^{-1}\left({\bf X}^{\intercal}{\boldsymbol{\psi}}\right).$ (14)

From here, one computes:

$\displaystyle e_{x}=\frac{12(1-\nu_{xy}\nu_{yx})d_{x}}{L_{z}^{3}},\qquad
e_{y}=\frac{12(1-\nu_{xy}\nu_{yx})d_{y}}{L_{z}^{3}},\qquad
g_{xy}=\frac{3d_{s}}{L_{z}^{3}},$ (15)

approximating $E_{x}$, $E_{y}$, $G_{xy}$. Given their small variability across
specimens, it will be assumed that the Poisson ratios $\nu_{xy}$, $\nu_{yx}$
are known and obtained from tabulated values.

Note that the methodology relies on the knowledge of the modal coefficients
$a_{m,n},b_{m,n},c_{m,n}$. These are, generally, unknown, but they can be
extracted by fitting the numerical eigenfrequencies of a number of training
plates, as shown in the next section.

##  3 Numerical Benchmark Tests

The methodology described above is now tested in a few numerical case studies.
First, the validity of the linear relationship (9) is assessed; second, the
ability of the current formulation to predict the correct rigidity constant
values from a batch of numerical test plates is shown.

###  3.1 Linear dependence of the non-dimensional frequencies on $p,q$

A key assumption in the current formulation is the linear dependence of the
non-dimensional frequencies $\bar{\omega}^{2}_{m,n}$ on the elastic constant
ratios $p,q$, as per (9). This is now checked by computing the numerical
eigenfrequencies of a set of plates created by varying reference plate
parameters. All throughout, the aspect ratio and the boundary conditions are
kept fixed: $\sigma=223/114$, $\mathcal{B}=\\{\text{C-F-F-F}\\}$, with the
clamped edge lying along $y$ (the radial direction). The finite element
simulation software package COMSOL is used for the numerical simulations. The
_Plate_ COMSOL component is selected. Note that the input COMSOL parameters
are dimensional and include the thick-plate parameters in addition to the
eight parameters listed in Section 2.1, though the former have a small
influence on the first few eigenmodes for thin and moderately thick plates.
Table 1 summarises the dimensional parameters used in the simulations.

Kirchhoff plate parameters
---
$\rho$ | $L_{x}$ | $L_{y}$ | $L_{z}$ | $E^{0}_{x}$ | $E^{0}_{y}$ | $G_{xy}$ | $\nu_{xy}$
473.9 | 22.3 | 11.4 | 3 | 10.7 | 0.716 | 0.50 | 0.51
kg $\cdot$ m-3 | cm | cm | mm | GPa | GPa | GPa |

Thick-plate parameters
---
$E_{z}$ | $G_{yz}$ | $G_{xz}$ | $\nu_{yz}$ | $\nu_{xz}$
0.39 | 0.023 | 0.62 | 0.45 | 0.50
GPa | GPa | GPa |  |

Table 1: Input constants for the COMSOL simulations, inspired from the table
in [11] (page 96). These constants are used to create a “training” set from
which the modal coefficients $a_{m,n}$, $b_{m,n},c_{m,n}$ are estimated. The
plates included in the training set have material and geometrical constants as
per the table, and where $(E_{x},E_{y})$ are selected as the twenty-five
elements of the set
$\mathcal{T}=\\{E_{x},E_{y}\,|\,E_{x}=(0.8,0.9,1.0,1.1,1.2)E_{x}^{0},\,E_{y}=(0.8,0.9,1.0,1.1,1.2)E_{y}^{0}\\}$.

From the table, twenty-five combinations of the elastic constant ratios $p,q$
are created, with $4.3<p<6.6$, $0.28<q<0.45$. The modal frequencies computed
in COMSOL are then sorted according to mode number and belong to the modes
(1,0), (1,1), (2,0), (1,2), (2,1), (2,2). In order to assign the correct modal
frequency to each modal shape when batch exporting from COMSOL, the modal
assurance criterion (MAC) can be used. Let ${\bf U},{\bf U}^{\prime}$ denote
the column vectors containing the sampled modal displacements of two modes.
The MAC index, a scalar value, is defined as [33]:

$\text{MAC}({\bf U},{\bf U}^{\prime}):=\frac{|{\bf U}^{\intercal}{\bf
U}^{\prime}|^{2}}{\|{\bf U}\|^{2}\,\|{\bf U}^{\prime}\|^{2}}.$ (16)

Values close to one indicate that the modal shapes ${\bf U},{\bf U}^{\prime}$
are highly correlated. In COMSOL, the sampled modal shapes can be batch-
exported and MAC indices can be computed against a reference set of
eigenshapes.

Figure 2: Modal Assurance Criterion (MAC) for the plate $\sigma=223/114$,
$\mathcal{B}=\\{\text{C-F-F-F}\\}$, under various choices of the parameters
$p,q$. A white square denotes a high correlation between the reference mode
shapes (denoted by R), and the batch of modal shapes for the indicated $p,q$,
ordered according to increasing frequency from 1 to 6. Modal crossings may
take place in some cases, such as in the second square $(p=4.36,q=0.40)$. Some
ambiguity exists in identifying modes (2,0) and (1,2), as seen in the last two
matrices. These cases are discarded when performing the fit. The reference set
of modes is selected as the batch $p=4.36,q=0.29$.

Figure 2 reports a few MAC matrices. Note that modal crossings may sometimes
take place. Furthermore, in many cases out of the training batch, the modes
(2,0) and (1,2) cannot be easily recognised, as one can see from the same
figure. This happens because the nodal lines are bent compared to the
reference $(x,y)$ axes: these modes deform considerably when changing $p,q$,
not allowing a clear identification of the modal indices $(m,n)$. This is a
property of the quarter-sawn plates reported in [20, p.84]. The cases for
which a clear modal identification is difficult are therefore discarded. See
also Figure 3.

Figure 3: Modal shapes for a plate with $\sigma=223/114$,
$\mathcal{B}=\\{\text{C-F-F-F}\\}$, $p=6.53$, $q=0.29$. Modes (1,2) and (2,0)
are marked with a star as assigning unambiguous modal indices may be difficult
in some cases. See also Figure 2.

The coefficients $a_{m,n}$, $b_{m,n}$, $c_{m,n}$ can now be estimated in
Matlab via a constrained linear fit, such that $a_{m,n}$, $b_{m,n}$, $c_{m,n}$
are non-negative, as per (9). The results are summarised in Table 2. Note that
the coefficients of determination $(R^{2})$ are very close to one in all
cases, highlighting the ability of the linear relationship (9) to reproduce
the correct modal frequency as a function of $p,q$. A visual representation of
the fit results for the modes is given in Figure 4.

| (1,0) | (1,1) | (2,0) | (1,2) | (2,1) | (2,2)
---|---|---|---|---|---|---
$a_{m,n}$ | 3.21 | 3.99 | 117.0 | 14.4 | 124 | 134
$b_{m,n}$ | 0.00 | 0.94 | 41.9 | 1720 | 15.6 | 1910
$c_{m,n}$ | 0.00 | 40.5 | 13.6 | 155 | 328 | 1260
$R^{2}$ | 1.00 | 1.00 | 0.999 | 0.999 | 1.00 | 1.00
Table 2: Linear fit results for a plate with $\sigma=223/114$, and
$\mathcal{B}=\\{\text{C-F-F-F}\\}$ (a cantilever plate clamped along the $y$
(radial) direction). The coefficients of determination $R^{2}$ are also
listed.  Figure 4: Linear dependence of the non-dimensional frequencies
$\bar{\omega}^{2}_{m,n}$ on the elastic constant ratios $p,q$. Here, the plate
has $\sigma=223/114$, and $\mathcal{B}=\\{\text{C-F-F-F}\\}$ (a cantilever
plate clamped along the radial direction). The dots represent the non-
dimensional modal frequencies as a function of $p,q$, and the planes are
plotted according to (9) using the values of $a_{m,n}$, $b_{m,n}$, $c_{m,n}$
listed in Table 2 .

###  3.2 Assessment of the method in the case $\mathcal{B}=\\{\text{C-F-
F-F}\\}$

The modal coefficients from Table 2 can be used to estimate the rigidity
constants of a batch of test numerical plates sharing the same boundary
conditions $\mathcal{B}=\\{\text{C-F-F-F}\\}$ and aspect ratio
$\sigma=223/114$. It is helpful to test the methodology for plates presenting
$p,q$ outside of the training ranges to assess the predicting power of the
fit. Hence, three test plates are considered, all sharing the same parameters
as per Table 1, but selecting $E_{x}$, $E_{y}$, so as to place $p,q$ inside or
outside their respective training ranges. Table 3 collects all such
parameters.

| $E_{x}/E_{x}^{0}$ | $E_{y}/E_{y}^{0}$ | $p$ | $q$ | $\text{err}_{x}(\%)$ | $\text{err}_{y}(\%)$ | $\text{err}_{s}(\%)$
---|---|---|---|---|---|---|---
Plate 1 | 1.3 | 0.7 | 7.1 | 0.25 | 0.77 | 0.9 | 0.78
Plate 2 | 1.0 | 0.7 | 5.4 | 0.25 | 0.65 | 1.2 | 0.031
Plate 3 | 1.3 | 1.0 | 7.1 | 0.36 | 0.28 | 0.11 | 0.57
Table 3: Test plate parameters and errors for the case $\sigma=223/114$,
$\mathcal{B}=\\{\text{C-F-F-F}\\}$. The other plate parameters are as per
Table 1. The training ranges, from Section 3.1, are: $4.3<p<6.6$,
$0.28<q<0.45$. The errors are defined as:
$\text{err}_{x}=100(D_{x}-d_{x})/D_{x}$,
$\text{err}_{y}=100(D_{y}-d_{y})/D_{y}$,
$\text{err}_{s}=100(D_{s}-d_{s})/D_{s}$. The modes considered in this test are
(1,0),(1,1),(2,1),(2,2). The values of $E^{0}_{x}$, $E^{0}_{y}$ are given in
Table 1.

For all test plates, the eigenfrequencies are computed numerically, and the
rigidity constants $d_{x}$, $d_{y}$, $d_{s}$ are estimated according to (14),
using the values of the modal coefficients ${\bf a},{\bf b},{\bf c}$ as in
Table 2.

The results are reported in Table 3, displaying errors below 1.2%. It is worth
commenting on the results. An important aspect concerns the selection of an
appropriate batch of modes for the estimation of the elastic constants. Modes
presenting an influence on all three constants may be used, provided that at
least three of them present linearly independent modal coefficients
$a_{m,n},b_{m,n},c_{m,n}$, as is the case for the coefficients listed in Table
2. From the same table, it is clear that using only modes (1,0), (1,1) and
(2,1) will not yield reliable results, as $D_{y}$ has no influence on these
modes’ dynamics (the $b_{n,m}$ coefficients are either zero or small compared
to $a_{n,m}$, $c_{n,m}$). In other words, the matrix in the least-square
formula (14) becomes poorly conditioned [34].

####  3.2.1 Sensitivity analysis for the case $\mathcal{B}=\\{\text{C-F-
F-F}\\}$

The numerical benchmark test carried out above highlights the ability of the
proposed method to recover the correct elasticity constants with very small
errors. However, the “experimental” frequencies were obtained via controlled
numerical simulation and are unaffected by all sources of error besides a
small inaccuracy due to the numerical approximation. Furthermore, the dynamics
of real plates may be somewhat more complicated than what is implied by model
(1): perfectly orthotropic conditions are seldom encountered in wooden boards;
thickness profiles may not be completely uniform; the boundary conditions,
especially of fixed type, may not be implemented exactly. Whilst care must be
taken to ensure that the experimental plate and the measurement setup are as
close to the ideal case as possible, experimental frequencies are likely
affected by some form of error. Thus, it is important to test the reliability
of the method using perturbed frequencies.

| $\epsilon_{(1,0)}$ | $\epsilon_{(1,1)}$ | $\epsilon_{(2,0)}$ | $\epsilon_{(1,2)}$ | $\epsilon_{(2,1)}$ | $\epsilon_{(2,2)}$ | $\text{err}_{x}(\%)$ | $\text{err}_{y}(\%)$ | $\text{err}_{s}(\%)$
---|---|---|---|---|---|---|---|---|---
Test 1 | 1.8 | -1.7 | 1.1 | -0.73 | 1.6 | 0.95 | -4.2 | 13.0 | -3.8
Test 2 | -1.6 | 0.76 | 1.6 | -1.1 | -1.4 | -2.0 | -3.5 | 7.2 | 11.0
Test 3 | 0.03 | -0.1 | 1.6 | -0.95 | 1.3 | -1.4 | -6.4 | 11.0 | 8.6
Test 4 | 0.76 | 0.89 | 0.66 | 0.79 | 1.7 | 1.6 | -3.4 | 9.2 | -6.0
Test 5 | -0.94 | -0.55 | 0.026 | -0.35 | 0.15 | -0.92 | -2.4 | 9.1 | 4.3
Table 4: Sensitivity test for the Plate 1 of Section 3.2, Table 3, with
$\sigma=223/114$, $\mathcal{B}=\\{\text{C-F-F-F}\\}$. The perturbation
$\epsilon$ was selected randomly within the interval $\pm 2\%$ with respect to
the exact modal eigenfrequency, and is reported in the table as a percentage.
The definitions of $\text{err}_{x}$, $\text{err}_{y}$ and $\text{err}_{s}$ are
as per Table 3.

Results are reported in Table 4, where the rigidity constants of Plate 1 from
Table 3 are computed using a perturbed set of eigenfrequencies, with
perturbations generated randomly within $\pm 2\%$ of the original
eigenfrequencies. As expected, the errors are higher, though still in line
with error bounds reported with previous methods, such as via the FEMU-3DVF
[3].

###  3.3 Assessment of the method in the case $\mathcal{B}=\\{\text{C-C-
C-C}\\}$

As a second case, consider the fully-clamped case with $\sigma=150/103$.
First, the modal coefficients $a_{m,n}$, $b_{m,n}$, $c_{m,n}$ are computed as
explained in Section 3.1.

| (0,0) | (0,1) | (1,0) | (0,2) | (1,1) | (1,2)
---|---|---|---|---|---|---
$a_{m,n}$ | 236 | 237.5 | 1790 | 242 | 1792 | 1799
$b_{m,n}$ | 1214 | 8589 | 1652 | 31910 | 10180 | 35250
$c_{m,n}$ | 146.9 | 552.3 | 534.9 | 1196 | 2051 | 4469
$R^{2}$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
Table 5: Linear fit results for a plate with $\sigma=150/103$, and
$\mathcal{B}=\\{\text{C-C-C-C}\\}$. The coefficients of determination $R^{2}$
are also listed.  Figure 5: Modal shapes for a plate with $\sigma=150/103$,
$\mathcal{B}=\\{\text{C-C-C-C}\\}$, $p=6.53$, $q=0.29$.

These are reported in Table 5. The coefficients of determination ($R^{2}$) are
all equal to 1.00 in this case, highlighting the validity of the linearity
assumption (9) in the fully-clamped case. The modal shapes for the same modes
are reported in Figure 5.

| $E_{x}/E_{x}^{0}$ | $E_{y}/E_{y}^{0}$ | $p$ | $q$ | $\text{err}_{x}(\%)$ | $\text{err}_{y}(\%)$ | $\text{err}_{s}(\%)$
---|---|---|---|---|---|---|---
Plate 1 | 1.3 | 0.7 | 7.1 | 0.25 | 0.81 | 0.86 | 0.12
Plate 2 | 1.0 | 0.7 | 5.4 | 0.25 | -0.28 | -0.26 | 0.043
Plate 3 | 1.3 | 1.0 | 7.1 | 0.36 | \- 0.12 | -0.12 | 0.036
Table 6: Test plate parameters and errors for the case $\sigma=150/103$,
$\mathcal{B}=\\{\text{C-C-C-C}\\}$. The training ranges, from Section 3.1,
are: $4.3<p<6.6$, $0.28<q<0.45$. The errors are defined as:
$\text{err}_{x}=100(D_{x}-d_{x})/D_{x}$,
$\text{err}_{y}=100(D_{y}-d_{y})/D_{y}$,
$\text{err}_{s}=100(D_{s}-d_{s})/D_{s}$. All six modes are considered in this
test.

A numerical benchmark test, similar to what was done for the previous case, is
reported in Table 6. The elastic constants are estimated with errors below
0.9%. A sensitivity analysis, analogous to what was done previously, is
performed on the fully clamped plate. The results, reported in Table 7, show
that $D_{x}$, $D_{y}$ are estimated consistently well across all tests, and
within small error bounds. A somewhat larger error is observed in the
estimation of $D_{s}$. In order to account for such a larger sensitivity, a
statistical approach can be adopted, as will be shown in Section 5.

| $\epsilon_{(0,0)}$ | $\epsilon_{(0,1)}$ | $\epsilon_{(1,0)}$ | $\epsilon_{(0,2)}$ | $\epsilon_{(1,1)}$ | $\epsilon_{(1,2)}$ | $\text{err}_{x}(\%)$ | $\text{err}_{y}(\%)$ | $\text{err}_{s}(\%)$
---|---|---|---|---|---|---|---|---|---
Test 1 | 0.64 | 0.51 | 0.37 | -1.0 | 1.6 | -0.95 | -2.5 | 3.7 | 10.6
Test 2 | 2.0 | -1.0 | -1.7 | -0.86 | -1.4 | 1.7 | 8.3 | 8.4 | -52
Test 3 | 1.2 | 1.4 | 2.0 | -1.9 | -0.65 | 0.63 | -1.1 | 5.7 | -6.9
Test 4 | -0.49 | -1.2 | -1.1 | 0.8 | 0.98 | 1.3 | 3.7 | 3.3 | -30.0
Test 5 | 0.89 | -0.95 | -0.2 | -1.9 | 0.75 | -1.6 | -1.1 | 5.4 | 11.0
Table 7: Sensitivity test for the Plate 1 of Table 6, with $\sigma=150/103$,
$\mathcal{B}=\\{\text{C-C-C-C}\\}$. The perturbation $\epsilon$ was selected
randomly within the interval $\pm 2\%$ with respect to the exact modal
eigenfrequency, and is reported in the table as a percentage. The definitions
of $\text{err}_{x}$, $\text{err}_{y}$ and $\text{err}_{s}$ are as per Table 6.

###  3.4 Interpretation of the results in terms of the Moore–Penrose inverse

In the classic methods by McIntyre and Woodhouse and by Caldersmith and
Freeman [1, 2], it is possible to extract the elasticity constants easily
since these appear decoupled in two “beam” modes and one “shear” mode of the
fully-free plate, $\mathcal{B}=\\{\text{F-F-F-F}\\}$. No such triplet of modes
exists for different boundary conditions, particularly of fixed type where the
effects of the elastic constants are densely coupled for all the modes, as one
can see easily from Table 5. Yet, the proposed methodology is able to retrieve
the correct elastic constant values with very small errors. The linear
formulation of the method offers a viable interpretation of such results. From
(14), one has that the rigidity constants are obtained in terms of the Moore-
Penrose inverse ${\bf X}^{\dagger}$, given by [35]:

${\bf X}^{\dagger}:=({\bf X}^{\intercal}{\bf X})^{-1}{\bf X}^{\intercal}.$
(17)

This linear operation reduces the overdetermined system comprising $N>3$
equations to a $3\times 3$ system, obtained as a linear combination of the $N$
available modes with the shortest Euclidian distance between the experimental
and reconstructed eigenfrequencies. It is natural to combine the eigenshapes
via ${\bf X}^{\dagger}$ and to observe the resulting modal shapes. This
operation is unnecessary to determine the elastic constants as such, though it
helps develop an intuitive interpretation of the results. The operation is a
simple matrix multiplication:

$\begin{bmatrix}{\boldsymbol{\alpha}}^{\intercal}\\\
{\boldsymbol{\beta}}^{\intercal}\\\
{\boldsymbol{\gamma}}^{\intercal}\end{bmatrix}={\bf
X}^{\dagger}\begin{bmatrix}{\bf U}_{1}^{\intercal}\\\ ...\\\ {\bf
U}_{N}^{\intercal}\end{bmatrix},$ (18)

where ${\boldsymbol{\alpha}}$, ${\boldsymbol{\beta}}$, ${\boldsymbol{\gamma}}$
are the combined mode shape column vectors, and where ${\bf U}_{i}$ is the
$i^{th}$ normalised column eigenshape. A visual representation is offered in
Figure 6, where the effects of $E_{x}$, $E_{y}$ and $G_{xy}$ are now decoupled
for modes $\alpha,\beta,\gamma$, and are clearly visible.

Figure 6: Combined mode shapes for the fully-clamped plate. These are obtained
by combining the mode shapes of Figure 5 using the Moore–Penrose inverse ${\bf
X}^{\dagger}$, so to obtain the modes $\alpha,\beta,\gamma$ for which the
dependence on, respectively, $E_{x}$, $E_{y}$ and $G_{xy}$ is completely
decoupled.

The possibility of combining modes in the current framework is one of the most
interesting findings of this work allowing the implementation of the method
using various choices of boundary conditions. The next section presents two
such experiments.

##  4 Experimental Setup

In order to validate the presented methodology, two experimental plates are
now considered, with the same boundary conditions and aspect ratios as the
numerical plates presented in Section 3. They are representative of
soundboards used in the production of stringed musical instruments: a kantele
and an acoustic guitar.

###  4.1 Tonewood specimen with $\mathcal{B}=\\{\text{C-F-F-F}\\}$

A sample of Finnish tonewood spruce was obtained from a master luthier based
in Finland. The sample was a leftover coming from the production of a concert
kantele, a Finnish plucked stringed musical instrument [36]. This board has
dimensions $L_{x}^{0}=25.3$ cm, $L_{y}=11.4$ cm, $L^{0}_{z}=5$ mm. The
thickness is reduced to $L_{z}=3$ mm using a belt sander, improving the thin-
plate approximation.

To realise the experimental boundary conditions, a clamping system is
implemented as in a previous study [37]. Two angular 1.8 kg iron elements are
clamped together by using a total of six spring plastic clamps distributed
along the entire length of the two angular components. To better distribute
the pressure along the $y$ axis, an additional flat 1.2 kg iron plate is
placed between one of the two angular elements and the specimen under test
while an additional layer of rubber material is wrapped around the clamped
plate’s side to prevent unwanted plate rotations. Figure 7 shows the final
experimental setup for the plate boundary conditions.

Figure 7: Experimental setup for the tonewood sample with
$\mathcal{B}=\\{\text{C-F-F-F}\\}$, $\sigma=223/114$, for which the modal
coefficients are listed in Table 2. The hammer is mounted on a pendulum which
stops after recoiling to avoid double hits.

The clamping mechanism shortens the specimen along the $x$ direction, such
that the effective length results as $L_{x}=22.3$ cm. The final plate
dimensions $L_{x}$, $L_{y}$, $L_{z}$ as well as the density $\rho$, are the
same as the ones reported in Table 1. In particular, the modal coefficients
$a_{m.n}$, $b_{m,n}$, $c_{m.n}$ for this plate are listed in Table 2.

In this experiment, a miniature impact hammer (PCB 086E80) is used to excite
the plate and a monoaxial accelerometer (PCB 352A21) is used to measure the
plate’s impulse response. A vinyl cover is used on the hammer’s hard tip to
reduce the damage caused by the impacts on the wooden plate surface. This
softer tip limits the force spectrum of the excitation signal [38], though it
still delivers sufficient energy to excite the board in the frequency range of
interest (below 1 kHz). To conduct the measurements, the impact hammer roves
across six equally spaced measurement positions whereas the accelerometer is
mounted and kept still in the same position during the entire measurement
campaign. To account for the measurement uncertainty caused by the low
reproducibility of the impact hammer force signal [39], five different impulse
responses are averaged over each measurement position. Particular care is
taken to avoid placing the accelerometer close to the nodal lines of the
interested modes listed in Table 2. Finally, the collected data is post-
processed in Matlab to analyse the frequency response functions from each
measurement position.

Figure 8 shows the resulting frequency spectra and Table 8 summarises the
experimental eigenfrequencies identified through a “peak finding” routine
implemented in Matlab and averaged across all the measurement positions.

Figure 8: Frequency spectra obtained from the six measurement points assessed on the tonewood spruce plate. The frequency range of interest here is up to 800 Hz. Dashed lines mark the spatially averaged experimental frequencies. | (1,0) | (1,1) | (2,0) | (1,2) | (2,1) | (2,2)
---|---|---|---|---|---|---
Measured freqs. (Hz) | 52 | 98 | 310 | 337 | 398 | 637
Table 8: Spatially averaged experimental frequencies and identified mode
shapes for the cantilever Finnish spruce tonewood.

As mentioned in Section 2, a crucial step of the presented methodology
requires the identification of the nodal lines along the $x$ and $y$
directions of the tonewood’s fibres, in order to associate the measured
frequencies with the corresponding modal coefficients
$a_{m,n},b_{m,n},c_{m,n}$. The modal shapes can be simply observed through the
well-known technique by Chladni [26]. Here, a mini shaker (PCB 2004E) is used
to excite the cantilever plate which is now horizontally oriented to
facilitate the mounting of the shaker to the plate surface (as shown in Figure
9) through the use of adhesive material (i.e. beeswax). Pure tones having the
same frequencies as the six measured eigenfrequencies reported in Table 8 are
used as excitation signals for the plate covered with off-the-shelf white
sand. Figure 9 shows the observed Chladni patterns.

Figure 9: Experimental setup including the miniature shaker (left) and Chladni
patterns (right) on the tonewood sample.

###  4.2 Acoustic guitar top board with $\mathcal{B}=\\{\text{C-C-C-C}\\}$

In a second experiment, a full guitar board is considered. This is the
assembly of two rectangular quarter-sawn spruce billets, referred to as a
_book-matched top_. Glueing together the half-boards to form a full guitar
board represents a preliminary step in the construction and design of guitars
conducted by luthiers and builders [40]. The board was purchased from
Ciresa111https://www.ciresafiemme.it/, a Fiemme Valley red spruce reseller,
and is here tested under fully clamped boundary conditions, mimicking those
occurring in the final assembled instrument [15]. To do so, an 8 kg plexiglass
frame is built and used to clamp the board along all its edges. Here, equally
distributed clamps are mounted along the frame’s edges and are used to apply
pressure on an additional plexiglass counter-frame which is resting directly
on the board under test. The experimental setup is shown in Figure 10. The
clamped board’s dimensions are: $L_{x}=60$ cm, $L_{y}=41.2$ cm, $L_{z}=3$ mm,
and its density is $\rho=399$ kg $\cdot$ m-3. Thus, the aspect ratio is
$\sigma=105/103$, and the modal coefficients $a_{m,n},b_{m,n},c_{m,n}$ are as
per Table 5.

Figure 10: Experimental setup and Chladni patterns for the Fiemme Valley red
spruce guitar top.

Given the larger size and weight of the guitar board compared to the
cantilever plate illustrated in Section 4.1, a higher force amplitude is
needed to excite all the interested modes of vibration with sufficient energy.
Furthermore, as the experimental fully clamped constraints require the board
to be now horizontally fixed to a frame, a different experimental setup is
here employed for ease of operation [41, 42]. Accordingly, the guitar board is
excited by an 8 $\Omega$ nominal impedance 48 g electro-dynamic exciter with a
25 mm diameter voice coil. Such a setup allows for a higher amplitude of
excitation and ease of operation on the horizontal board.

An exponential sine sweep ranging from 45 Hz to 8 kHz is used as the test
signal. Similarly to the previous measurement campaign, the exciter roves
across four different measurement points distributed over the guitar’s board.
As a receiver, an Earthworks MD30 class-1 microphone is placed at a 125 mm
distance from the board to capture the output signal from the guitar board.
Finally, the recorded data is post-processed using the Aurora plugins [43],
yielding the impulse responses which are then further processed and analysed
in Matlab. Table 9 reports the measured eigenfrequencies.

| (0,0) | (0,1) | (1,0) | (1,1) | (0,2) | (1,2)
---|---|---|---|---|---|---
Measured freqs. (Hz) | 57 | 100 | 130 | 162 | 170 | 224
Table 9: Spatially averaged experimental frequencies and identified mode
shapes for the fully clamped red spruce plate.

For the presented case study, Chladni patterns are obtained by exciting the
plate at the measured frequencies using pure tones and dried tea leaves, and
are visible in Figure 10.

##  5 Results and Discussion

The elastic constants for the Finnish and the Fiemme Valley spruce tonewoods
are retrieved using the experimental frequencies from Tables 8 and 9
respectively, the modal coefficients from Tables 2 and 5, and formulae (14)
and (15). Various combinations of modes can be used here. Both cases include
six identified experimental modes; thus, it is possible to run (14) using
forty-two combinations containing at least three modes out of six. This allows
performing many estimates of the elastic constants, yielding a statistically
significant set from which mean values and standard deviations can be
computed. It is best to exclude from the statistics all spurious estimates.
First, negative estimates should be excluded (these appear in some cases when
the Moore-Penrose is poorly conditioned [34], as discussed previously).
Second, outliers should be excluded. These may be given in terms of deviation
from the mean: here, all estimates lying in the twentieth percentile away from
the mean are excluded. Results are summarised in Figures 11 and 12.

Figure 11: Elastic constant estimates for the Finnish spruce tonewood. Forty-
two mode sets are considered, corresponding to all possible combinations
containing at least three experimental modes out of the six measured ones.
Mean values and standard deviations are represented as dashed lines and
coloured bands, respectively. Outliers are selected as belonging to the
twentieth percentile away from the mean. Figure 12: Elastic constant estimates
for the Fiemme Valley red spruce tonewood. The same analysis as per Figure 11
holds.

The mean values and standard deviations of $e_{x}$, $e_{y}$, $g_{xy}$ for the
two plates are summarised in Table 10. Note that the standard deviations are
very low in all cases (below 5%), except for the shear modulus in the fully
clamped case. This is a consequence of the sensitivity of the method under
fully-clamped conditions, as discussed in Section 3.3, though the statistical
approach adopted here is able to contain the deviation to within a reasonable
threshold. The estimated mean values are in line with previously reported
values for spruce, such as in [11, p.96].

|  |  $e_{x}$ (GPa) |  $e_{y}$ (MPa) |  $g_{xy}$ (MPa)
---|---|---|---|---
Finnish spruce (cantilever) | mean | 13.0 | 702 | 668
std | 2.8% | 2.2% | 4.1%
Fiemme Valley spruce (clamped) | mean | 10.4 | 1003 | 530
std | 1.73% | 3.7% | 14%
Table 10: Elastic constants mean values and standard deviations for the two
spruce tonewoods. Values are in line with previously reported values for
spruce, such as in [11, p.96].

The average elastic constant values can be used as input parameters in COMSOL,
and the numerical frequencies obtained can be assessed against the measured
experimental frequencies. This is done in Table 11 and 12, where the numerical
frequencies are in both cases predicted very accurately: errors fall below
2$\%$ in all cases, highlighting the accuracy of the proposed methodology in
retrieving appropriate elastic constant values.

| Meas. (Hz) | Num. (Hz) |  $\Delta f_{n}$ (Hz) |  $\frac{\Delta f_{n}}{f_{n}}$ (%) |  ${\Delta f_{n}}$ (cent)
---|---|---|---|---|---
(1,0) | 52 | 51 | -1 | -1.5 | -25
(1,1) | 98 | 100 | 2 | 1.9 | 32
(2,0) | 311 | 315 | 4 | 1.3 | 22
(1,2) | 337 | 336 | -1 | -0.4 | -7
(2,1) | 398 | 394 | -4 | -0.9 | -16
(2,2) | 637 | 629 | -8 | -1.2 | -21
Table 11: Errors between the experimentally measured and numerically computed eigenfrequencies, for the Finnish spruce plate, under cantilever boundary conditions. The numerical frequencies were computed in COMSOL using the mean elastic constant values from Table 10.  | Meas. (Hz) | Num. (Hz) |  $\Delta f_{n}$ (Hz) |  $\frac{\Delta f_{n}}{f_{n}}$ (%) |  ${\Delta f_{n}}$ (cent)
---|---|---|---|---|---
(0,0) | 57 | 56 | -1 | -1.8 | -31
(0,1) | 100 | 98 | -2 | -2 | -35
(1,0) | 130 | 130 | 0 | 0 | 0
(1,1) | 162 | 161 | -1 | -0.6 | -11
(0,2) | 170 | 169 | -1 | -0.6 | -10
(1,2) | 224 | 221 | -3 | -1.3 | -23
Table 12: Errors between the experimentally measured and numerically computed
eigenfrequencies, for the Fiemme Valley red spruce plate, under fully clamped
boundary conditions. The numerical frequencies were computed in COMSOL using
the mean elastic constant values from Table 10.

##  6 Conclusion

In this work, a method to estimate the elastic constants of thin, rectangular
orthotropic panels was offered, as an application of non-destructive inverse
parameter estimation methods. Compared to previously available techniques in
the literature, the proposed framework does not rely on a specific set of
boundary conditions or modes, allowing multiple estimates of the elastic
constants via simple matrix inverses. This way, not only can mean values be
estimated, but also deviations. The enabling idea is the identification of a
linear relationship between the elastic constants and the square of the non-
dimensional plate frequencies, with modal-dependent coefficients that can be
easily computed from a batch of numerical training plates. An interpretation
of the method was offered in terms of the Moore-Penrose inverse, which is used
to combine the eigenmodes to yield three modes in which the influence of the
thin-plate elastic constants is completely decoupled. A number of numerical
and experimental tests showed the reliability of the proposed methodology in
cases of interest in musical acoustics.

## Acknowledgements

This work received funding from the European Research Council (ERC) under the
Horizon2020 framework, grant number 950084 - StG - NEMUS. Francesco Vai and
Roberto Budini from the Department of Industrial Engineering at the University
of Bologna are kindly acknowledged for designing and building the impact
hammer pendulum. Ciresa is kindly acknowledged for assisting in the selection
of the Fiemme Valley tonewood. Master luthier Pekka Lovikka is thanked for
providing the Finnish tonewood samples.

## References

  * [1] M.E. McIntyre and J. Woodhouse.  On measuring the elastic and damping constants of orthotropic sheet materials.  Acta Metall, 36(6):1397–1416, 1988.
  * [2] G. Caldersmith and E. Freeman.  Wood properties from sample plate measurements. Part I.  J Catgut Acoust, 1(2):8–12, 1990.
  * [3] R. Viala, V. Placet, and S. Cogan.  Identification of the anisotropic elastic and damping properties of complex shape composite parts using an inverse method based on finite element model updating and 3D velocity fields measurements (FEMU-3DVF): Application to bio-based composite violin soundboards.  Compos Part A-Appl S, 106:91–103, 2018.
  * [4] K. Ege, N. B. Roozen, Q. Leclere, and R. G. Rinaldi.  Assessment of the apparent bending stiffness and damping of multilayer plates; modelling and experiment.  J Sound Vib, 426:129–149, 2018.
  * [5] S. Jin, D. Gruber, and H. Harmuth.  Determination of Young’s modulus, fracture energy and tensile strength of refractories by inverse estimation of a wedge splitting procedure.  Eng Fract Mech, 116:228–236, 2014.
  * [6] S. Jin, D. Gruber, H. Harmuth, and R. Rössler.  Thermomechanical failure modeling and investigation into lining optimization for a Ruhrstahl Heraeus snorkel.  Eng Fail Analysis, 62:254–262, 2016.
  * [7] A. Maciolek, R. Wagener, and T. Melz.  Review of and a new approach to elastic modulus evaluation for fatigue design of metallic components.  Int J Fatigue, 151:106325, 2021.
  * [8] Z. Chen, U. Gandhi, J. Lee, and R.H. Wagoner.  Variation and consistency of Young’s modulus in steel.  J Mat Process Tech, 227:227–243, 2016.
  * [9] J. H. Tam, Z. C. Ong, Z. Ismail, B. C. Ang, and S. Y. Khoo.  Identification of material properties of composite materials using nondestructive vibrational evaluation approaches: A review.  Mech Adv Mat Struct, 24(12):971–986, 2017.
  * [10] E. Brühwiler and F.H. Wittmann.  The wedge splitting test, a new method of performing stable fracture mechanics tests.  Eng Fract Mech, 35(1):117–125, 1990.
  * [11] V. Bucur.  Handbook of materials for string musical instruments.  Springer, Cham, Switzerland, 2016.
  * [12] J. E. Mottershead, M. Link, and M. I. Friswell.  The sensitivity method in finite element model updating: A tutorial.  Mech Syst Signal Proc, 25(7):2275–2296, 2011.
  * [13] A. van Casteren, W.I. Sellers, S.K.S. Thorpe, S. Coward, R.H. Crompton, and A.R. Ennos.  Factors affecting the compliance and sway properties of tree branches used by the Sumatran orangutan (Pongo Abelii).  PLOS ONE, 8(7):e67877, 2013.
  * [14] K. R. Bootle.  Wood in Australia : types, properties and uses.  McGraw-Hill, 2010.
  * [15] T. Gore and G. Gilet.  Contemporary Acoustic Guitar Design and Build: Build.  Trevor Gore, 2011.
  * [16] R. You, N. Zhu, X. Deng, J. Wang, and F. Liu.  Variation in wood physical properties and effects of climate for different geographic sources of Chinese fir in subtropical area of China.  Sci Rep-UK, 11(1):4664, 2021.
  * [17] M. Albano, D. Comelli, G. Fiocco, M. Mattonai, J. J. Lucejko, L. Zoia, M. P. Colombini, and M. Malagodi.  Chemical modification of wood induced by the traditional making procedures of bowed string musical instruments: the effect of alkaline treatments.  Herit Sc, 10(1):76, 2022.
  * [18] R. Malvermi, M. Albano, S. Gonzalez, G. Fiocco, F. Antonacci, M. Malagodi, and A. Sarti.  The impact of alkaline treatments on elasticity in spruce tonewood.  Sci Rep-UK, 12(1):13335, 2022.
  * [19] P. Persson and O. Flodén.  Effect of material parameter variability on vibroacoustic response in wood floors.  Appl Acoust, 146:38–49, 2019.
  * [20] N. H. Fletcher and T. D Rossing.  The physics of musical instruments.  Springer Science & Business Media, 2012.
  * [21] Voichita Bucur.  Acoustics of wood.  Springer Science & Business Media, 2006.
  * [22] R. Malvermi.  A statistical approach to non-destructive testing in stringed musical instruments.  PhD thesis, Politecnico di Milano, 2023.
  * [23] M. Quintavalla, F. Gabrielli, and C. Canevari.  Grading materials for stringed instruments soundboards: An approach considering the orthotropic elastic and damping properties.  Appl Acoust, 187:108521, 2022.
  * [24] J. V. Burke and M. C. Ferris.  A Gauss—Newton method for convex composite optimization.  Math Progr, 71(2):179–194, 1995.
  * [25] Thomas D Rossing.  Chladni’s law for vibrating plates.  Am J Phys, 50(3):271–274, 1982.
  * [26] Ernst Florens Friedrich Chladni.  Treatise on Acoustics: The First Comprehensive English Translation of EFF Chladni’s Traité d’Acoustique.  Springer, 2015.
  * [27] G. Aksu.  Dynamic analysis of orthotropic plates using a finite difference formulation.  PhD thesis, Loughborough University, 1974.
  * [28] A. Wilczyński and M. Kociszewski.  Determination of elastic constants of particleboard layers by compressing glued layer specimens.  Wood Res, 56(1):77–91, 2011.
  * [29] S.S. How, H.S. Sik, and I. Ahmad.  Review on six types of log cutting methods in various applications: Part I.  Timber Tech Bull, 45, 2007.
  * [30] S. Berg, E. Unsal, and H. Dijk.  Non-uniqueness and uncertainty quantification of relative permeability measurements by inverse modelling.  Comput Geotech, 132:103964, 2021.
  * [31] F. Esqueda, B. Kuznetsov, and J. D. Parker.  Differentiable white-box virtual analog modeling.  In 24th International Conference on Digital Audio Effects (DAFx), pages 41–48, Vienna, Austria, September 2021.
  * [32] S. P Boyd and L. Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * [33] M. Pastor, M. Binda, and T. Harčarik.  Modal assurance criterion.  Procedia Eng, 48:543–548, 2012.
  * [34] Y. Wei, H. Diao, and S. Qiao.  Condition number for weighted linear least squares problem.  J Comput Math, pages 561–572, 2007.
  * [35] J. C. A. Barata and M. S. Hussein.  The Moore–Penrose pseudoinverse: A tutorial review of the theory.  Brazilian J Phys, 42:146–165, 2012.
  * [36] H. Tahvanainen.  On the acoustics of the concert kantele.  In Stockholm Music Acoustics Conference, Stockholm, Sweden, June 2023.
  * [37] S. Duran, M. Ducceschi, H. Tahvanainen, and L. Ausiello.  Experimentally-tuned synthesis of a thin plate.  In Proceedings of the Institute of Acoustics, Winchester, UK, October 2023.
  * [38] D.J. Ewins.  Modal testing: theory, practice and application.  John Wiley & Sons, 2009.
  * [39] P. Avitabile.  Modal testing: A practitioner’s guide.  John Wiley & Sons, 2017.
  * [40] M. Rau.  Measurements and analysis of acoustic guitars during various stages of their construction.  J Acoust Soc Am, 149(4 Supplement):A25–A25, 2021.
  * [41] L. Ausiello, L. Yule, G. Squicciarini, and C. Barlow.  Guitar soundboard measurements for repeatable acoustic performance manufacturing.  In Reproduced Sound, Bristol, UK, November 2018. Institute of Acoustics.
  * [42] L. Ausiello and V. Hockey.  Quantitative measurements to enhance performance of acoustic musical instruments and improve manufacturing.  Acoustics Bulletin, 47(2), 2021.
  * [43] A. Farina.  Advancements in impulse response measurements by sine sweeps.  In Audio Engineering Society convention 122, Vienna, Austria, May 2007. Audio Engineering Society.