# Non parametric harmonization of a tabla's membrane

Sylvain Maugeais

► To cite this version:

Sylvain Maugeais. Non parametric harmonization of a tabla's membrane. 2023. hal-04102160

HAL Id: hal-04102160

https://hal.science/hal-04102160

Preprint submitted on 22 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Non parametric harmonization of a tabla's membrane

Sylvain Maugeais

## Abstract

The musical drums of the Indian subcontinent are known for their harmonic quality due to their non uniform membrane, requiring a precise distribution of mass. Extending the works of previous authors, the present work uses a non-parametric optimisation procedure to find an optimal density for a membrane loaded by air. The modes of the loaded membrane are computed using a Finite Element Method coupled with a Boundary Element Method, and a minimum of inharmonicity if found using a gradient algorithm with an explicitly computed gradient.

## 1 Introduction

The tabla, along with other Indian percussion instruments like pakhawaj and mridangam, is known amongst membranophones for its harmonic quality (cf. [18]) and resonance (cf. [19], chapter 6.3). These specificities are mainly due to the syahi, which is an extra mass added at the center of the membrane. Although it used to be made of clay that had to be applied each time the instrument was used (cf. [3] chapter 33, verse 25-26), it is now made either of single use dough (bass side of pakhawaj of mridangam) or syahi: "a paste of starch, gum, iron oxide, charcoal, or other materials." (cf. [8] chapter 18.5) (for tablas). The process of either applying the dough (cf. [3]) or the syahi (cf. [19]) is quite complex. The latter requires the application of different layers of syahi paste, drying after each step and rubbing with a stone. These different steps ensure there is a reticulum of cracks forming in the syahi (see fig. (1)), which reduces the stiffness of the added mass (cf. [19], chapter 5.11).

<figure><img src="image_2.png" /><figcaption>Figure 1: Close up of the syahi. The reticulum of cracks and different layers can be seen</figcaption></figure>

The shape of the syahi has to be very precise to achieve harmonicity. Leaving aside the sustain, various authors have tried to describe the its geometry by numerically optimising the inharmonicity: from [17] (refined by [10]) viewing it as a composite membrane, [21] and [20] in a similar way but with a smooth approximation of the composite density, to [16] and [1] who used non paramet-

ric optimisation. All these articles found inhomegeneous membranes whose first in vacuo eigenfrequencies are almost harmonic, and whose shape is very similar to that of an actual physical tabla.

On the other hand, following [5] who studied air loading on timpani membrane, [22] applied the same optimisation procedure as [17] on composite membranes but taking into account the air loading. They too obtained an almost harmonic set of frequencies. As a byproduct, they also got decay times. Unfortunately, these are sometimes unrealistic mainly because their model does not take into account the viscoelastic losses, which can be incorporated along the line of [9] who studied it for the kettledrum.

The present article expands these results by applying a non parametric optimisation (like [1]), taking into account air loading (like [22]) and viscoelastic losses (like [9]), and imposing conditions on the density ensuring the playability of the tabla (like [16]).

Following [12] in their work on optimisation of musical instruments, the finite element method has been chosen to compute the eigenfrequencies relative to general densities. Together with boundary element method and optimisation, these techniques form the technical core of the present work.

The remainder of the article is organised as follows: in section 2, the model used for the tabla and air loading is described and the modes in vacuo are computed using the finite element method (FEM). In section 3, the matrices involved in the boundary element method (BEM) are computed, first outside the kettle, and then inside. In section 4, the solutions of the complete system, membrane together with air loading, is computed, and are compared with the results of the bibliography. Finally, in section 5, the actual optimisation is carried out and discussed. The article concludes with open questions and future directions.

## 2 Problem formulation

Throughout the text, the only concern is that of the right hand tabla. The left hand one, although it has also been studied in other works (eg. [22]), is much more difficult to model as the player's wrist is usually used during playing by resting at different positions and applying different pressures while fingers strike the membrane ([13] chapter
<table><caption>Table 1: List of symbols</caption>
  <tr>
    <td>Σ</td>
    <td>Surface of the drumhead</td>
  </tr>
  <tr>
    <td>F</td>
    <td>Infinite flange</td>
  </tr>
  <tr>
    <td>a</td>
    <td>Radius of the drumhead (in m) = 0.05m for a tabla</td>
  </tr>
  <tr>
    <td>L</td>
    <td>Height of the kettledrum (in m) = 0.1m for a tabla</td>
  </tr>
  <tr>
    <td>x = (ρ, θ, z)</td>
    <td>Point in space, cylindrical coordinates</td>
  </tr>
  <tr>
    <td>r = (ρ, θ)</td>
    <td>Point in the membrane, polar coordinates</td>
  </tr>
  <tr>
    <td>σ</td>
    <td>Surface density of the drumhead (in kg/m²) = 0.245kg/m² for a uniform tabla membrane</td>
  </tr>
  <tr>
    <td>η</td>
    <td>Transverse displacement of the drumhead (in m)</td>
  </tr>
  <tr>
    <td>T</td>
    <td>Membrane tension (in N/m) = 1822N/m for a tabla</td>
  </tr>
  <tr>
    <td>ν</td>
    <td>Viscoelastic coefficient of the drumhead (in s) = 0.6 × 10⁻⁶s</td>
  </tr>
  <tr>
    <td>p<sub>in</sub></td>
    <td>Pressure inside the kettledrum</td>
  </tr>
  <tr>
    <td>p<sub>out</sub></td>
    <td>Pressure above the flange F</td>
  </tr>
  <tr>
    <td>ρ<sub>0</sub></td>
    <td>Air density = 1.20kg/m³</td>
  </tr>
  <tr>
    <td>c<sub>a</sub></td>
    <td>Speed of sound = 344m/s</td>
  </tr>
</table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Model of the tabla</figcaption></figure>

7), therefore changing considerably the boundary conditions, and producing a continuous range of pitches.

## 2.1 Model

Following [5], the tabla is modeled as a rigid cylinder of radius $a$ and height $L$, capped at the top by a membrane $Σ$ and at the bottom by a rigid disk. An infinite rigid plane flange enclosing the membrane in the plane of the membrane is assumed (see ibid. concerning the validity of this assumptions).

Let us write $\mathbf{x} = (r, \theta, z)$ for the cylindrical coordinates, the membrane $Σ$ being positioned in the $z = 0$ plane and the center of the membrane at $r = 0$.

Following [5] (and [9] for the viscoelastic term), the displacement $\eta$ of the membrane is governed by the equation

$$ \sigma \frac{\partial^2 \eta}{\partial t^2} = T \Delta \left( \eta + \nu \frac{\partial \eta}{\partial t} \right) + p_{in}|_{z=0} - p_{out}|_{z=0} \quad (1) $$

where $\sigma(\mathbf{x})$ is the (non uniform) density of the membrane at the point $\mathbf{x}$, $T$ is the membrane tension by unit length, $\nu$ is the viscoelastic damping coefficient, $p_{in}$ and $p_{out}$ are the acoustic pressure fields inside and outside the kettle (i.e. above the plane $z = 0$). Assuming the membrane is fixed at the rim yields

$$ \eta(\mathbf{x}, t) = 0 \text{ for all } t \text{ and } \mathbf{x} \in \partial\Sigma \quad (2) $$

As the kettle acts as a baffle, a Neuman condition is used on the fixed borders of the kettle

$$ \left. \frac{\partial p_{in}}{\partial \rho} \right|_{\rho=a} = \left. \frac{\partial p_{in}}{\partial z} \right|_{z=-L} = 0 \text{ for all } t. \quad (3) $$

In the same way, the assumption that the membrane is surrounded by an infinite plane flange implies

$$ \left. \frac{\partial p_{out}}{\partial \rho} \right|_{\rho>a, z=0} = 0 \quad \text{for all } t. \quad (4) $$

Denote by $\hat{\eta}$ (resp. $\hat{p}_{in}, \hat{p}_{out}$) the temporal Fourier transform of $\eta$ (resp. $p_{in}, p_{out}$) with variable $\omega$. The global acoustic pressure field $p$ composed of $p_{in}$ and $p_{out}$ satisfies the wave equation

$$ \left( \Delta + \frac{\omega^2}{c_a^2} \right) \hat{p} = 0 \quad (5) $$

where $c_a$ denotes the speed of sound in air.

Denote by $G_{in}^{\omega}$ and $G_{out}^{\omega}$ the Green functions for equation (5) for the inside and outside domains. Then the pressure fields can be computed (cf. [5], eqs. (21) and (22))

$$ \hat{p}_{in}(\mathbf{x}) = \int_{\Sigma} G_{in}^{\omega}(\mathbf{x}|\mathbf{x}') \left. \frac{\partial \hat{p}_{in}}{\partial z'} \right|_{z'=0} d\mathbf{x}' \quad (6) $$

$$ \hat{p}_{out}(\mathbf{x}) = \int_{\Sigma} G_{out}^{\omega}(\mathbf{x}|\mathbf{x}') \left. \frac{\partial \hat{p}_{out}}{\partial z'} \right|_{z'=0} d\mathbf{x}'. \quad (7) $$

The compatibility of the speeds of air and membrane at $z = 0$, together with linearized fluid dynamics, imply

$$ \left. \frac{\partial \hat{p}_{in}}{\partial z} \right|_{z=0} = \left. \frac{\partial \hat{p}_{out}}{\partial z} \right|_{z=0} = \omega^2 \rho_0 \hat{\eta}, \quad (8) $$

with $\rho_0$ denoting the density of air at equilibrium. Finally the equation verified by $\hat{\eta}$ is given by

$$ - \omega^2 \sigma \hat{\eta} = T \Delta (\hat{\eta} - j \omega \nu \hat{\eta}) + \omega^2 \rho_0 \left( \int_{\Sigma} G_{in}^{\omega}(\mathbf{x}|\mathbf{x}') \left. \hat{\eta}(\mathbf{x}') \right|_{z'=0} d\mathbf{x}' - \int_{\Sigma} G_{out}^{\omega}(\mathbf{x}|\mathbf{x}') \left. \hat{\eta}(\mathbf{x}') \right|_{z'=0} d\mathbf{x}' \right). \quad (9) $$

## 2.2 Description of the elements and solutions in vacuo

The membrane is meshed uniformly with triangles and a spectral element method of order $n$ is used. Let $S$ be the set of shape functions defined on the triangulation using Legendre-Gauss-Lobatto points that are not on the boundary of $\partial\Sigma$. Define a hermitian product on the Sobolev space $H_0^1(\Sigma)$ by

$$ \langle f, g \rangle = \int_{\Sigma} f(\mathbf{x}) \overline{g(\mathbf{x})} d\mathbf{x}. \quad (10) $$

For the density function $\sigma$ on $\Sigma$ and tension parameter $T$, define the matrices

$$ \mathbf{A}_{\sigma} = (\langle \sigma \varphi, \psi \rangle)_{\varphi, \psi \in S} \quad (11) $$

$$ \mathbf{B} = T (\langle \operatorname{grad} \varphi, \operatorname{grad} \psi \rangle)_{\varphi, \psi \in S}. \quad (12) $$

The computation of the scalar products for computing the matrices $\mathbf{A}_{\sigma}$ and $\mathbf{B}$ are performed using the collocation method.
<figure><img src="image_4.png" /><figcaption>Figure 3: Relative error in cents for the modes in vacuo for a uniform circular membrane compared to theoretical values, for different uniform meshes and different orders.</figcaption></figure>

The lossless in vacuo equation implied by (9) is given by

$$-\omega^2 \mathbf{A}_{\sigma} \Phi = \mathbf{B} \Phi \quad (13)$$

and its solutions $\omega, \Phi$ are obtained through the computation of eigen-pairs, and are Galerkin weak approximations of in vacuo frequencies and modes of the membrane with density $\sigma$.

According to [2] (see also [15] for a general discussion on the size of the mesh), the error term for $\omega$ or $\Phi$ can be bounded for small enough eigenvalues by $h^n\omega^2$, where $h$ denotes the mesh-size, which is the maximum diameter of the triangles of the mesh.

The $\omega$s are counted with multiplicities and sorted by their magnitude as a set $(\omega_i)_{i \in I}$ where $I$ denotes the set of eigen-pairs of (13), and the corresponding eigenvector is denoted by $\Phi_i$.

The error on frequencies in cent for different meshes and different orders are given on figure (3) for the first 120 modes (extra logarithmic scale on the cents is used to allow the comparison between high and low order). This number of modes is necessary to ensure a good precision during the next step of the computation (cf. section 4.1). As expected, the error increases with frequency (cf. [2]) and should be minimal for the first ones. The poor precision on the first eigenvalues is explained by the approximation of the geometry: a lesser number of points in the mesh implies a worse approximation of the circle.

As the number of modes is important for the precision required in section 4.1, and the speed of computation is not relevant, order 3 and 1096 points are used in the rest of the article.

The numerical method was implemented in Python from scratch on a desktop computer without using any finite element package.

# 3 Boundary element method

Contrary to [5], [22] and [9], the integrals in (9) are not computed “explicitly” because the density of the membrane is not constant, so the projection onto the eigenfunctions of the radial component of the membrane equation cannot be readily applied. Instead, the Boundary Element Method (BEM) is used (cf. for example [11] or [14]), together with a spectral element methods (SEM) of order 3. All the integrals involving only regular functions are computed using collocation methods.

## 3.1 Computation of the outside Green function

According to [5], equation (20), $G_{out}$ can be expressed as

$$G_{out}^{\omega}(\mathbf{x}|\mathbf{x}') = \frac{-1}{4\pi} \left( \frac{e^{j \frac{\omega}{c_a} \|\mathbf{x}-\mathbf{x}'\|}}{\|\mathbf{x}-\mathbf{x}'\|} + \frac{e^{j \frac{\omega}{c_a} \|\mathbf{x}^*-\mathbf{x}'\|}}{\|\mathbf{x}^*-\mathbf{x}'\|} \right) \quad (14)$$

with $\mathbf{x}^* = (r, \theta, -z)$ if $\mathbf{x} = (r, \theta, z)$. In particular, when $\mathbf{r}, \mathbf{r}' \in \Sigma$ we get

$$G_{out}^{\omega}(\mathbf{r}|\mathbf{r}') = \frac{-2}{4\pi} \frac{e^{j \frac{\omega}{c_a} \|\mathbf{r}-\mathbf{r}'\|}}{\|\mathbf{r}-\mathbf{r}'\|} \quad (15)$$

Define the matrix $\mathbf{G}_{out}(\omega)$ whose entries are

$$(\mathbf{G}_{out}^{\omega})_{\varphi, \psi} = \frac{-1}{4\pi} \int_{\Sigma} \varphi(\mathbf{r}) \int_{\Sigma} \psi(\mathbf{r}') \frac{e^{j \frac{\omega}{c_a} \|\mathbf{r}-\mathbf{r}'\|}}{\|\mathbf{r}-\mathbf{r}'\|} d\mathbf{r}' d\mathbf{r}. \quad (16)$$

for any two $\varphi, \psi \in S$ shape functions.

Denote by $\kappa_{\psi}(\mathbf{r}) = \int_{\Sigma} \psi(\mathbf{r}') \frac{e^{j \frac{\omega}{c_a} \|\mathbf{r}-\mathbf{r}'\|}}{\|\mathbf{r}-\mathbf{r}'\|} d\mathbf{r}'$ so that

$$(\mathbf{G}_{out}^{\omega})_{\varphi, \psi} = \frac{-1}{2\pi} \int_{\Sigma} \varphi(\mathbf{r}) \kappa_{\psi}(\mathbf{r}) d\mathbf{r} \quad (17)$$

For each $s \in S$, denote by $\mathbf{r}_s$ its node point and $w_s$ its weight in the Lobatto quadrature. As $\kappa_{\varphi}$ is regular, $(\mathbf{G}_{out}^{\omega})_{\varphi, \psi}$ is approximated by

$$(\mathbf{G}_{out}^{\omega})_{\varphi, \psi} \approx \frac{-1}{4\pi} \sum_{s \in S} w_s \varphi(\mathbf{r}_s) \kappa_{\psi}(\mathbf{r}_s) = w_{\varphi} \kappa_{\psi}(\mathbf{r}_{\varphi}) \quad (18)$$

To compute $\kappa_{\psi}$, the usual technique is used (cf. [6], 6.1): separate the singular part and compute it analytically, the rest of the integral is computed numerically. For this write $\frac{e^{j \frac{\omega}{c_a} \|\mathbf{r}_{\varphi}-\mathbf{r}'\|}}{\|\mathbf{r}_{\varphi}-\mathbf{r}'\|} = \frac{1}{\|\mathbf{r}_{\varphi}-\mathbf{r}'\|} - E(\|\mathbf{r}_{\varphi}-\mathbf{r}'\|, \omega)$, with $E$ is a regular function. Then

$$\kappa_{\psi}(P_{\varphi}) = \underbrace{\int_{\Sigma} \psi(\mathbf{r}') \frac{1}{\|\mathbf{r}_{\varphi} - \mathbf{r}'\|} d\mathbf{r}'}_{\kappa_{\varphi, \psi}^*} + \int_{\Sigma} \psi(\mathbf{r}') E(\|\mathbf{r}_{\varphi} - \mathbf{r}'\|, \omega) d\mathbf{r}'$$

The first integral $\kappa_{\varphi, \psi}^*$ can be computed explicitly using a polar parametrization of triangles, together with antiderivatives of trigonometric polynomials. As for the
<figure><img src="image_5.png" /><figcaption>Figure 4: Relative error in cents for the eigenvalues of the Laplace operator with zero Neumann boundary condition, for different uniform meshes and different orders.</figcaption></figure>

second, E being regular, it can be approximated using the collocation method

$$\int_{\Sigma} \psi(\mathbf{r}') E_{\varphi}(\mathbf{r}', \omega) d\mathbf{r}' \cong w_{\psi} E(\|\mathbf{r}_{\varphi} - \mathbf{r}_{\psi}\|, \omega).$$

Finally

$$(\mathbf{G}_{\text{out}}^{\omega})_{\varphi, \psi} \cong \frac{-1}{2\pi} w_{\varphi} (\kappa_{\varphi, \psi}^{*} + w_{\psi} E(\|\mathbf{r}_{\varphi} - \mathbf{r}_{\psi}\|, \omega))$$

which is an approximation of order n thanks to Lobatto's quadrature.

Finally, define the matrices **H** and **E**($\omega$) by

$$\mathbf{H}_{\varphi, \psi} = w_{\varphi} \kappa_{\varphi, \psi}^{*} \quad (19)$$

$$\mathbf{E}_{\varphi, \psi}(\omega) = w_{\varphi} w_{\psi} E(\|\mathbf{r}_{\varphi} - \mathbf{r}_{\psi}\|, \omega) \quad (20)$$

so that

$$\mathbf{G}_{\text{out}}^{\omega} = \frac{-1}{2\pi}(\mathbf{H} + \mathbf{E}(\omega)) \quad (21)$$

## 3.2 Computation of the Green function inside the shell

Let $(\Psi_n)_{n \in \mathbb{N}}$ be a complete system of eigenvectors for the Laplace operator on $\Sigma$ with Neumann boundary conditions, and denote by $\varpi_n$ its eigenvalue. This data can be computed using the finite element method as above changing the boundary condition, the weak approximation now taking place in the Sobolev space $H^1(\Sigma)$ instead of $H_0^1(\Sigma)$.

Figure 4 gives the error in the computations of the $\varpi_n$ for different meshes and different orders. As for the Dirichlet condition, order 3 with 1096 points give very good results.

A computation following [5] equation (19) gives an explicit form for the Green function inside the shell:

$$G_{in}^{\omega}(\mathbf{r}|\mathbf{r}') = - \sum_{n} \Psi_{n}(\mathbf{r}) \Psi_{n}(\mathbf{r}') \frac{\cotan(\mu_{n}L)}{\mu_{n}} \quad (22)$$

with $\mu_n = \sqrt{\varpi_n^2 - \frac{\omega^2}{c_a^2}}$.

To compute the weak approximation matrix of $G_{in}^{\omega}$ at entries $\varphi, \psi \in S$ remark that

$$(\mathbf{G}_{in}^{\omega})_{\varphi, \psi} = \int_{\Sigma} \int_{\Sigma} G_{in}^{\omega}(\mathbf{r}|\mathbf{r}') \varphi(\mathbf{r}) \psi(\mathbf{r}') d\mathbf{r} d\mathbf{r}', \quad (23)$$

so that

$$(\mathbf{G}_{in}^{\omega})_{\varphi, \psi} = - \sum_{n} \frac{\cotan(\mu_n L)}{\mu_n} \langle \Psi_n, \varphi \rangle \langle \Psi_n, \psi \rangle \quad (24)$$

## 3.3 Projection onto the space of in vacuo modes

To keep track of the right symmetries when looking for eigenmodes, as well as reducing the dimension of the search space, it is convenient to use a projection of the equations onto the space of eigenmodes in vacuo $H_{vac} = \text{vect}(\Phi_i)_{0 \le i \le l}$. It is a good candidate for the approximation of the exact solution u thanks to Sturm-Liouville's theory. Restricting to this space is the equivalent of the truncation of the infinite series done in [5] for equations (38) and (46).

For any matrix **N** among **A**, **B**, $\mathbf{G}_{in}^{\omega}$, $\mathbf{G}_{out}^{\omega}$, define its restriction to the space $H_{vac}$ by

$$\tilde{\mathbf{N}}_{i,j} = {}^t \bar{\Phi}_i \mathbf{N} \Phi_j. \quad (25)$$

for all $0 \le i, j \le l$.

The projection of the fundamental equation (9) can then be written as

$$-\omega^2 \tilde{\mathbf{A}}_{\sigma} M = (1 - j\omega\nu) \tilde{\mathbf{B}} M + \omega^2 \rho_0 (\tilde{\mathbf{G}}_{in}^{\omega} - \tilde{\mathbf{G}}_{out}^{\omega}) M \quad (26)$$

In general, the results depend on the value of $l$, a bigger value of $l$ giving more precise results. We found that a typical value of $l = 80$ is a good choice for the computations to ensure a good precision for the first 15 modes (cf. section 5). This high value of $l$ is the main reason why order 3 methods were chosen (see figures 3 and 4).

# 4 Solutions of the complete system

The last step to find the eigenmodes of the loaded non uniform membrane is to solve equation (26) which is non polynomial in $\omega$ because of the presence of this variable in the matrices $\tilde{\mathbf{G}}_{in}^{\omega}$ and $\tilde{\mathbf{G}}_{out}^{\omega}$. This is achieved following the same “fixed point” strategy as in [5]. It should be noted however that this method does not always give a solution (see ibid. section IV) and is also confronted to the well known problem of spurious eigenvalues and non-uniqueness of solutions (see [4]). Although this is a real problem for bigger drums, we found that all the solutions given by this method for the tabla are of the same order as that of a physical tabla, owing to the relatively small size of the membrane and the reduced air loading compared to that of the timpani.
## 4.1 Fixed point algorithm

The method chosen here follows only the “simple case” of [5] section IV: start with an eigensolution in vacuo ($\omega_0, M_0$), where $M_0$ denotes the coordinates of the solution in the basis of the $\Phi_t$s. In particular, $M_0$ is part of the canonical basis:

$$M_0 = {}^t(0, \dots, 0, \overset{m}{1}, 0, \dots, 0).$$

Denote by $m$ the non zero coordinate.

Then solve recursively the eigenvalue problem

$$-\zeta^2 \tilde{\mathbf{A}}_{\sigma} \Lambda = \underbrace{(1 - j\omega_n \nu) \tilde{\mathbf{B}} + 2\omega_n^2 \rho_0 (\tilde{\mathbf{G}}_{in}^{\omega_n} - \tilde{\mathbf{G}}_{out}^{\omega_n})}_{\text{depends only on } \omega_n} \Lambda \quad (27)$$

with indeterminates $\zeta$ and $\Lambda$, where the previous guess $\omega_n$ has been replaced in the Green functions and the viscoelastic losses.

In general, this system has $\ell + 1$ eigensolutions. To ensure that the chosen solution has the same kind of symmetries as the initial guess $M_0$ in terms of nodal circles, the one that has the largest relative $m$-th coordinate (hence closest to $M_0$) is picked. This is the crucial step and requires that the perturbation is not too strong ([5] section IV). This hypothesis not only allows the fixed point method to converge, but also gives realistic results for the tabla.

The fixed point method is described precisely below (cf. algorithm 1).

<caption>**Algorithm 1:** Fixed point algorithm</caption>

**Data:** $m < \ell$
**Result:** solution of (26)
$\bar{\omega}_0 = \omega_m$

**do**
| $S \leftarrow \{(\zeta_s, M_s)\}$ = set of normalized eigen-solutions of (27)
| $i \leftarrow \operatorname*{argmax}\{m$-th coordinate of$|M_s|, (\zeta_s, M_s) \in S\}$
| $\bar{\omega}_{n+1} \leftarrow \zeta_i$
| $n \leftarrow n + 1$
**while** $|\bar{\omega}_{n+1} - \bar{\omega}_n| > \varepsilon$

The output of this algorithm is compared on figure 5 for the uniform membrane of a timpani to the results of [5] table I, for different meshes and different orders. The modes (0, 1) and (0, 2) where removed because the algorithm does not converge in that case (or sometimes converges to a spurious mode) which is in agreement with the discussion of [5], section IV. Indeed, these two modes are the most perturbed ones, and the real solutions are far from the modes in vacuo. However, for the other modes, the results are within the margin of error as the results in op. cit. are given with a precision of only 10 cents. The results obtained in this figure confirm that a high order scheme is necessary, together with a fine mesh.

<figure><img src="image_6.png" /><figcaption>Figure 5: Relative error in cents for the timpani compared to [5], table I, for different uniform meshes and different orders. Modes (0, 1) and (0, 2) where removed because of convergence problems.</figcaption></figure>

## 4.2 Comparison with known results

### Case of tabla

The results obtained by the fixed point algorithm for different uniform meshes and different orders are compared to those of [22], table II for the in vacuo modes (see figure 6) and for modes loaded by air (see figure 7). As for the uniform membrane detailed on figure 5, the results are within the margin of error, as the results in [22] Table II are given only with a precision of 10 cents. However, the effect of high order is less clear on these figures than for figure 5. This is due to the discontinuity of the composite membranes which induces less regularity on the FEM problem and a slower convergence.

The modes in vacuo and loaded are well approximated at order 2 and 3 for a uniform mesh of 1096 points.

### Regularisation of the composite membrane

To get a clearer result, the method is tested on a regularisation of the composite density and compared to the computation with the same algorithm but with degree 3 and finer mesh (1096 points). The results can be seen on figure 8 for different mesh sizes and are compatible with convergence of the method.

### On the necessity of the viscoelastic losses

Because the only dissipative mechanism used in the models of [5] and [22] is the sound radiation, both articles give unrealistic values for $\tau_{60}$'s of some modes (cf. [5], section IV and [22], section III.C). To get more realistic decay times, a straightforward improvement is to add viscoelastic losses, as in [9] who studied the timpani model of [5] together with such an added term. From here onwards, we consider the same value as [9] for the viscoelastic coefficient: $\nu = 0.6 \times 10^{-6}$s.
<figure><img src="image_7.png" /></figure>

Figure 6: Relative error in cents for inharmonicities of the modes in vacuo relative to mode (0, 2), compared to [22], table II, for different uniform meshes and different orders. Mode (0, 2) is removed as it is taken as reference

<figure><img src="image_8.png" /></figure>

Figure 8: Relative error in cents for the eigenfrequencies of the loaded modes for a regularisation of the composite density of [22], table II, for four different uniform meshes. The comparison is made with order 3 and 1096 points. The top right corner shows the regularisation used

<figure><img src="image_9.png" /></figure>

Figure 7: Relative error in cents for inharmonicities of the modes in vacuo relative to mode (0, 2) compared to [22], table II, for four different uniform meshes and different orders. Mode (0, 2) is removed as it is taken as reference

<figure><img src="image_10.png" /></figure>

Figure 9: Comparison between the $\tau_{60}$ obtained by [5], table III (black curve for the computation, blue curve for the experiment) and the FEM/BEM method (green curve).
The comparison of $\tau_{60}$'s obtained from the FEM/BEM method exposed in section 3 taking the viscoelastic losses into account are presented on figure 9. As pointed out in [9], the addition of the viscoelastic losses give more realistic values for all the modes and the remaining differences could come from the model of the room acoustics (cf. [22]).

# 5 Optimisation

The goal of this section is to give the optimisation tools necessary to find a solution to the tuning process. The method chosen here relies on the gradient algorithm as a “natural” way to go toward a minimum. This method is not only suitable because of the big dimension space in which the optimisation process takes place (the dimension of the space is given by the number of forms in the FEM method, which is here about 10 000) but also because it gives a natural, local direction to search for the minimum.

## 5.1 Computation of gradients

The goal of this section is to study the dependence of solutions of the system

$$ -\sigma\omega^2\hat{\eta} = (1 - j\nu\omega)T\Delta\hat{\eta} + G_{\omega} * \hat{\eta} \quad (28) $$

under small deformations of the functional parameter $\sigma$. For this, consider the perturbed system

$$ \begin{aligned} -(\sigma + \varepsilon\tilde{\sigma})(\omega + \varepsilon\tilde{\omega})^2(\hat{\eta} + \varepsilon\tilde{\eta}) = & \\ & (1 - j\nu(\omega + \varepsilon\tilde{\omega}))T\Delta(\hat{\eta} + \varepsilon\tilde{\eta}) + \\ & \rho_0(\omega + \varepsilon\tilde{\omega})^2 G_{\omega + \varepsilon\tilde{\omega}} * (\hat{\eta} + \varepsilon\tilde{\eta}) \end{aligned} \quad (29) $$

where $\varepsilon$ denotes a “small” parameter.

Taking the Taylor expansion in $\varepsilon$ and keeping only the part in $\varepsilon$ give

$$ \begin{aligned} & -(\tilde{\sigma}\omega^2\hat{\eta} + 2\sigma\omega\tilde{\omega}\hat{\eta} - j\nu\tilde{\omega}T\Delta\hat{\eta} \\ & + \rho_0\tilde{\omega}\left(2\omega G_{\omega} + \omega^2\frac{\partial G_{\omega}}{\partial\omega}\right)\hat{\eta}) = \\ & \qquad \sigma\omega^2\tilde{\eta} + T(1 - j\nu\omega)\Delta\tilde{\eta} + \rho_0\omega^2 G_{\omega}\tilde{\eta}. \end{aligned} \quad (30) $$

To simplify this equation, take the scalar product $\langle \cdot, \tilde{\eta} \rangle$ with the complex conjugate $\tilde{\eta}$ of $\hat{\eta}$ of the whole equation, and use the hermitian symmetries of the operators $\Delta$, $G_{\omega}$ and $\frac{\partial G_{\omega}}{\partial\omega}$. The fact that $(\omega, \hat{\eta})$ is a solution of (9) yields

$$ \begin{aligned} & \langle \tilde{\sigma}\omega^2\hat{\eta} + 2\sigma\omega\tilde{\omega}\hat{\eta} - j\nu\tilde{\omega}T\Delta\hat{\eta} + \\ & \qquad \rho_0\tilde{\omega}\left(2\omega G_{\omega} + \omega^2\frac{\partial G_{\omega}}{\partial\omega}\right) * \hat{\eta}, \tilde{\eta} \rangle = 0. \end{aligned} \quad (31) $$

Thereby giving

<figure><img src="image_11.png" /><figcaption>Figure 10: Evolution of frequencies along the line $\sigma_U + \delta\tilde{\sigma}$ together with the derivative computed with the real part of formula (32).</figcaption></figure>

$$ \text{grad}_{\sigma} \omega = \frac{-\omega^2}{2\omega \langle \sigma \hat{\eta}, \tilde{\eta} \rangle - j\nu T \langle \Delta \hat{\eta}, \tilde{\eta} \rangle + \rho_0 \omega \langle (2G_{\omega} + \omega \frac{\partial G_{\omega}}{\partial \omega}) * \hat{\eta}, \tilde{\eta} \rangle} \hat{\eta}^2 \quad (32) $$

Any perturbation of $\sigma$ along a direction $\varepsilon\tilde{\sigma}$ can then be computed as

$$ \omega(\sigma + \varepsilon\tilde{\sigma}) = \omega(\sigma) + \varepsilon\langle\tilde{\sigma}, \text{grad}_{\sigma}\omega\rangle \quad (33) $$

Note that this is true for both real and imaginary parts, so that it is possible to see the tendency of both frequency and decay time in a particular direction. This is what is done on figures 10 (for the frequency) and 11 (for decay time) taking for $\sigma$ the uniform density $\sigma_U = 0.245$, for $\tilde{\sigma}$ a regular function that cancels on the rim

$$ \tilde{\sigma}(\rho, \theta) = \sigma_U(1 - \rho/a)^2 $$

and a not so small parameter $\delta \in [-0.3, 0.3]$ so that the evolution along $\delta\tilde{\sigma}$ can be compared with the derivative computed using the gradient. Note that $\delta$ is also chosen so that the density remains positive. For both figures, the colored curves represent the frequency or decay time for the different modes computed with de FEM/BEM method (using a 1096 points mesh and order 3, together with $l = 80$) and the dashed lines represent the derivatives computed with formula of the gradient (formulas (32) and (33)).

Evolution of both frequencies and decay times are very well predicted by the gradient alone.

## 5.2 Definition of inharmonicity

As was already noticed by Raman in [18], the first nine modes of a tabla are in an almost harmonic progression. The different modes, defined by their nodes, are described on figure 12 together with the approximate ratio to the mode (1, 1) which serves as a reference.
<figure><img src="image_12.png" /><figcaption>Figure 11: Evolution of $\tau_{60}$ along the line $\sigma_U + \delta\tilde{\sigma}$ together with the derivative computed with the imaginary part of formula (32)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 13: Gradient of inharmonicity computed for a uniform density</figcaption></figure>

It should be noted that this description of the modes hides the fact that most of them are actually degenerate: for symmetry reasons, each mode $(\alpha, \beta)$ is actually doubled if $\alpha > 0$. However, among the twin modes, only the mode shapes differ (by a rotation) and not the frequency.

The inharmonicity function is therefore defined as

$$I(\sigma) = \sum_{m \in \text{modes}} \gamma_m \left( \frac{2\Re(\omega_m(\sigma))}{p_m \Re(\omega_{(1,1)}(\sigma))} - 1 \right)^2 \quad (34)$$

where $\gamma_m$ designates a weight affected to mode $m$. The choice of these weights is actually crucial as they favor which modes are to be optimised first.

For example, [21] chose $\gamma_m = p_m^2$ which gives a relatively strong weight to high frequency modes. This choice leads in our case to the gradient algorithm to not converge.

Another interesting choice is the one by [22], equation (24), which is closely related to the uniform case but taking into account psychoacoustic phenomena.

The present article follows [10] (with a different choice of target ratios $p_{(0,3)}$) and [16], giving the same weight to all modes $\gamma_m = 1$. This choice makes the gradient algorithm converge toward a physically realistic solution (cf. figure 17).

## 5.3 Monotonification

As can be seen on figure 13, the gradient of the inharmonicity is neither of constant sign, nor monotonic (as a function of $r$), and does not even have a precise circular symmetry. All these characteristics are however expected from the density of the syahi.

The non constancy of the sign is problematic to build the syahi, as it would necessitate to remove some matter, even maybe to the point of having negative density at some places: the syahi is usually built by adding small layers of mass.

The non circularity of the gradient comes from the fact that only a finite set of modes taken into account for the calculation of inharmonicity, together with the choice of the limits in the definition of inharmonicity and the precision of the computations. This “problem” is more of a numerical nature, but is relatively natural for the “monotonification” process described hereafter.

Finally, even for densities with circular symmetry, an optimal solution could have non playable shapes (for example those proposed in [1] Figure 7) which could be

<figure><img src="image_14.png" /><figcaption>Figure 12: First nine modes of the tabla and ratios $2f/f_{(1,1)}$</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: Function $\sigma_0 - \delta\text{grad}_\sigma I(\sigma_0)$ as a multivalued function of $r$ (in black) together with its monotonification $(\sigma_0 - \delta\text{grad}_\sigma I(\sigma_0))^{\sharp}$ (in red).</figcaption></figure>

harmful to the player. Moreover, owing to the current building process of the *syahi* (cf. [19] or [7]), which is done by adding small concentric layers of paste of decreasing radius, letting them dry and rubbing them with a stone to produce the characteristic cracks (cf. figure 1), one hypothesis that seems reasonable, and which is found in actual physical tablas, is to suppose that the density is monotonic as a function of *r*.

A simple way to achieve all three hypotheses is to use a process called here “monotonification” which is done in three steps: for any real valued function $\lambda$ on $\Sigma$

* transform $\lambda$ into a non negative function,

* take the maximum of $\lambda$ for every fixed radius $r$, thus ensuring the circular symmetry,

* take the cumulative maximum, thus ensuring the monotonicity in $r$.

These steps can be summarized in the formula

$$ \lambda^{\sharp}(r) = \max_{r' \ge r} \left( \max_{|\mathbf{r}|=r'} \left( \max(\lambda(\mathbf{r}), 0) \right) \right) \quad (35) $$

giving a new function which depends only on $r$, is non negative and decreasing.

On figure 14, the monotonified version of

$$ (\sigma_0 - \delta\text{grad}_\sigma I(\sigma_0))^{\sharp} $$

is given for a uniform density $\sigma = 0.245$ and a parameter $\delta$ given by the first step of the gradient algorithm (see algorithm 2). It should be noted that $\sigma_0 - \delta\text{grad}_\sigma I(\sigma_0)$ is drawn only as a function of $r$, which is why it appears as a multivalued function.

## 5.4 Gradient algorithm

The algorithm used in this article is the usual gradient descent algorithm. Although there are more effective algorithms to find a global minimum, this one is very suitable for our application as the dimension of the search space

is quite big (about 10 000) and it gives a result in only a few steps (about 30 steps for the order 3 FEM/BEM method with a mesh of 1096 points). The algorithm used for optimisation is detailed hereafter in algorithm 2.

<caption>**Algorithm 2:** Gradient algorithm with monotonification of $\sigma$</caption>

**Data:** $\sigma$
**Result:** locally optimal density
$n \leftarrow 0$
$\sigma_0 \leftarrow \sigma$
$I_0 \leftarrow \text{inharmonicity}(\sigma_0)$
**do**
| $\text{grad} = \text{gradient of inharmonicity}(\sigma_n)$
| $\delta \leftarrow -I_n / \langle \text{grad}, \text{grad} \rangle$
| $\sigma_{n+1} \leftarrow (\sigma_n + \delta \cdot \text{grad})^{\sharp}$
| $I_{n+1} \leftarrow \text{inharmonicity}(\sigma_{n+1})$
| $n \leftarrow n + 1$
**while** $I_n > I_{n-1}$

The only significative difference with the usual gradient algorithm is the use of monotonification in the step $\sigma_{n+1} \leftarrow (\sigma_n + \delta \cdot \text{grad})^{\sharp}$.

It should be noted that another choice is possible by applying the monotonification on the gradient, rather than on the density. This latter choice has the effect of producing a series of increasing densities: $\sigma_{n+1} \ge \sigma_n$, which may be closer to the reality of construction of the *syahi*, but give results that are further from the optimal.

## 5.5 Results of optimisation

The results of the gradient algorithm are given on figures 15 to 18 for an initial uniform density $\sigma_0 = 0.245 \text{kg/m}^2$ as in [22], Table II. The parameters used for the table are therefore $T = 1822N/m$, $a = 0.05m$ and $L = 0.1m$. Denoting by $\sigma_T$ the optimal density obtained in loc. cit., we get $I(\sigma_T) = 0.00229$ with formula 34, which is different from the inharmonicity computed in [22] because of the choice of weights. It should be noted that, although $\sigma_T$ was obtained for another choice of weights, it is very close to the optimal inharmonicity, emphasizing the robustness of the composite density.

The decrease of inharmonicity is clearly seen on figure 15 and the algorithm stops after 36 steps for an optimal density $\sigma_{opt}$ for which $I(\sigma_{opt}) = 0.00019$. The fact that $I(\sigma_T) > I(\sigma_{opt})$ is not a surprise as $\sigma_T$ verifies all the hypotheses (symmetry, monotonicity) that are imposed upon the search space, which therefore include that of [22].

The evolution of the frequency of each mode, normalised by the frequency of mode (1, 1), can be seen on figure 16. It seems that the higher the frequency, the slower the convergence. This is readily explained by the choice of weights that is made, which favours the low frequency.

The evolution of the density is displayed in figure 17, together with the composite optimal $\sigma_T$ in red for reference. The proximity of $\sigma_T$ and $\sigma_{opt}$ is striking, emphasizing the fact that seeing the *syahi* as a composite membrane is a very good approximation.
<figure><img src="image_16.png" /><figcaption>Figure 15: Evolution of inharmonicity during optimisation (blue curve). Inharmonicity of the optimal obtained in [22] is printed for reference (dotted red curve)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16: Evolution of inharmonicity of each mode relative to the mode (1, 1) during optimisation</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Evolution of density during optimisation, together with $\sigma_T$ in red for reference</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18: Evolution of $\tau_{60}$ of each mode during optimisation</figcaption></figure>

It should be noted that from the very beginning of the optimisation, the rightmost point where $\sigma_n > \sigma_0$ shifts to the left during optimisation. This is due to the choice of taking the monotonification on $\sigma$, rather than on $\text{grad}_{\sigma}I$. It allows the optimisation process to get a smaller support for the syahi.

A third ramp can be seen on the last obtained densities around $r = 0.01m$. This suggests that to improve even further the composite membrane of [22], another ring should be added (see [23] for a study in this direction without the air loading).

The evolution of decay time for each mode can be seen on figure 18. During all the optimisation process, the decay times are increasing, giving for the optimal density a sustain that is up to 10 times greater for $\sigma_{opt}$ than for $\sigma_0$, in particular for the low frequency modes. This is compatible with the knowledge of tabla makers, that adding smaller central circle adds to the resonance (cf. [19], section 6.3).
Sounds obtained at each step of the simulation
are available online at http://perso.univ-lemans.fr/~smauge/tabla/

## 6 Conclusion

A non parametric optimisation procedure is proposed to
produce a non homogeneous membrane loaded by air
whose first partials are harmonic. Together with mild hy-
potheses ensuring the playability of the membrane, this
procedure gives a realistic density that is close to that
of an actual physical tabla. As a byproduct, it gives an
explicit formula for the changes in frequencies and de-
cay times of eigenmodes with respect to an infinitesimal
variation of the density. Moreover, being based on finite
elements and boundary elements methods, it can readily
be applied to non circular geometries.

The method developed here gives a density close to
those obtained in previous studies, highlighting the good
aproximation given by a composite membrane. The main
novelty of the present article lies in the numerical scheme,
which uses a non parametric optimisation together with
mild hypotheses that are inherent to the manufacture of
the syahi : circular symmetry and radial monotonicity.
As such, it could be used in more complex settings where
no obvious solution such as the composite one is known.

As for the previously developed methods, the optimisa-
tion is done only for the first five frequencies (equivalent
to 15 eigenmodes, some of them being degenerate). This
is justified as these are the only prevalent frequencies ob-
served in the spectrum of the tabla. It was conjectured by
Raman in [18] that the reason for this limited spectrum is
the presence of a second, annular shaped membrane that
rests on the main one. We hope that the tools developed
in the present article will help investigate this conjecture.

One of the main drawback of the method is the assump-
tion that the stiffness is negligible. This would require a
precise investigation as the added mass is important and
the material appears to be quite rigid. As stated in the
introduction, the usual explanation for this low stiffness is
the reticulum of cracks which, to our knowledge, has not
been studied in depth and would deserve more attention.

References

[1] P. Antunes. Harmonic configurations of non-homogeneous membranes. *Acta Acustica united with Acustica*, 2005, 50, 2017.

[2] I. Babuška and J.E. Osborn. Eigenvalue problems. *In Handbook of Numerical Analysis, Vol. II*, page 641–787, 1991.

[3] Bharata-Muni. *The Natyasastra*. Bibliotheca Indica, Calcuta, 1951. Translated by Manmohan Ghosh.

[4] J.T. Chen. Recent development of dual bem in acoustic problems. *Computer Methods in Applied Mechanics and Engineering*, 188(4):833–845, 2000. IVth World Congress on Computational Mechanics. (II). Optimum.

[5] R.S. Christian, R.E. Davis, A. Tubis, C.A. Anderson, R.I. Mills, and T.D. Rossing. Effects of air loading on timpani membrane vibrations. *J. Acoust. Soc. Am.*, 76(5):1336–1345, 1984.

[6] M. Costabel. Principles of boundary element methods. *Computer Physics Reports*, 6(1):243–274, 1987.

[7] D. Courtney. *Manufacture and repair of tabla*. Sur Sangeet Services, Houston, 2001.

[8] N.H. Fletcher and Rossing D.T. *The physics of musical instruments*. Springer, Houston, 1998.

[9] E.U. Gallardo, M.A. Alonso-Arévalo, E. García-Canseco, and C. Aguilar-Ibáñez. Sound model of an orchestral kettledrum considering viscoelastic effects. *Applied Acoustics*, 164, 2020.

[10] S. Gaudet, C. Gauthier, and S. Léger. The evolution of harmonic indian musical drums: A mathematical perspective. *Journal of Sound and Vibration*, 291(1):388–394, 2006.

[11] W.S. Hall. *The Boundary Element Method*. Springer, Dordrecht, 1994.

[12] L. Henrique and J. Antunes. Optimal design and physical modelling of mallet percussion instruments. *Acta Acustica united with Acustica*, 89(6):948–963, 2003.

[13] J. Kippen. *The tabla of Lucknow : a cultural analysis of a musical tradition*. Cambridge studies in ethnomusicology. Cambridge University Press, Cambridge, 1988.

[14] S. Kirkup. The boundary element method in acoustics: A survey. *Applied Sciences*, 9(8):1642, 2019.

[15] S. Marburg. Six boundary elements per wavelength: is that enough? *Journal of computational acoustics*, 10(01):25–51, 2002.

[16] S. Maugeais. How to apply a plaster on a drum to make it harmonic. In *Proceedings of the International Symposium on Musical Acoustics*, 2014.

[17] B.S. Ramakrishna and M. M. Sondhi. Vibration of indian musical drums regarded as composite membranes. *J. Acoust. Soc. Am.*, (26):523–529, 1954.

[18] C.V. Raman. Indian musical drums. *Proc. Indian Acad. Sci., Sect. A*, (1A):179–188, 1934.

[19] P. A. Roda. *Resounding Objects: Musical Materialities and the Making of Banarasi Tablas*. PhD thesis, New York University, 2013.

[20] T. Samejima and R Fukuda. Vibration analysis of a musical drum head under nonuniform density and tension using a spectral method. *Acoustical Science and Technology*, 37:295–302, 2016.

[21] G. Sathej and R. Adikari. the eigenspectra of indian muscial drums. *J. Acoust. Soc. Am.*, (125 (2)):831–838, 2009.
[22] S. Tiwari and A. Gupta. Effects of air loading on the acoustics of an indian musical drum. *J. Acoust. Soc. Am.*, 141(4), 2017.

[23] G. Vautour, A. Brahmi, and C. Gauthier. Quasi-harmonic circular membranes with a central disc or an isolated ring of different density. *Journal of Sound and Vibration*, 332(19):4732-4740, 2013.