# Total harmonic distortion of a piezoelectric MEMS loudspeaker in an IEC 60318-4 coupler estimation using static measurements and a nonlinear state space model

Romain Liechti, Stéphane Durand, Thierry Hilt, Fabrice Casset, Christophe Poulain, Gwenaël Le Rhun, Franklin Pavageau, Hugo Kuentz, Mikaël Colin

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

## ► To cite this version:

Romain Liechti, Stéphane Durand, Thierry Hilt, Fabrice Casset, Christophe Poulain, et al.. Total harmonic distortion of a piezoelectric MEMS loudspeaker in an IEC 60318-4 coupler estimation using static measurements and a nonlinear state space model. Micromachines, 2021, 12 (12), pp.1437-10.3390/mi12121437. cea-04951765

HAL Id: cea-04951765

https://cea.hal.science/cea-04951765v1

Submitted on 17 Feb 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

# Total Harmonic Distortion of a Piezoelectric MEMS Loudspeaker in an IEC 60318-4 Coupler Estimation Using Static Measurements and a Nonlinear State Space Model

Romain Liechti<sup>1,2,*</sup>,<sup>1</sup> Stéphane Durand<sup>2,†</sup>, Thierry Hilt<sup>1</sup>, Fabrice Casset<sup>1</sup>, Christophe Poulain<sup>1</sup>,
Gwenaël Le Rhun<sup>1</sup>, Franklin Pavageau<sup>1</sup>, Hugo Kuentz<sup>1</sup> and Mikaël Colin<sup>1</sup>

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

check for
updates

**Citation:** Liechti, R.; Durand, S.; Hilt, T.; Casset F.; Poulain C.; Le Rhun G.; Pavageau F.; Kuentz H.; Colin M. Total Harmonic Distortion of a Piezoelectric MEMS Loudspeaker in an IEC 60318-4 Coupler Estimation Using Static Measurements and a Nonlinear State Space Model. *Micromachines* **2021**, *12*, 1437. https://doi.org/10.3390/mi12121437

Academic Editor: Libor Rufer

Received: 27 October 2021

Accepted: 20 November 2021

Published: 24 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** We propose a method to evaluate the Total Harmonic Distortion generated by a cantilever-based PZT loudspeaker inside an IEC 60318-4 coupler. The model is validated using experimental data of a commercial loudspeaker. Using the time domain equations of the equivalent electrical circuit of the loudspeaker inside the coupler and a state space formulation, the acoustic pressure response is calculated and compared to the measurement of the manufacturer. Next, the stiffness, transduction and capacitance nonlinear functions are evaluated with a Double-Beam Laser Interferometer (DBLI) and a nanoindenter on test devices and on the commercial loudspeaker. By introducing the nonlinear functions into the model as amplitude-dependent parameters, the THD generated by the loudspeaker is calculated and compared to the value provided by the manufacturer. The good agreement between the measurement and the simulation could allow for a rather quick simulation of the performance of similarly designed loudspeakers at the early stage of the design, by only estimating the static linearity of the main nonlinearity sources.

**Keywords:** MEMS; piezoelectric; loudspeaker; total harmonic distortion; state space; IEC 711; thin film; PZT

## 1. Introduction

Loudspeakers are used to converting an electrical signal into sound waves, as accurately as possible, with a sufficiently high sound pressure level. This is usually achieved with a piston-like membrane, behaving as a first-order oscillating system, actuated by an electromechanical transducer. Driven by the forever growing industry of portable and connected devices, research was carried out in order to make the loudspeaker compatible with microfabrication processes, in order to reduce the manufacturing tolerances, the cost and the thickness of the devices. An interesting performance was demonstrated with Micro Electromechanical Systems (MEMS) electrodynamic loudspeakers [1,2] and piezoelectric loudspeakers [3–7]. Nevertheless, the achieved performance is not competitive compared to non-MEMS loudspeaker. However, as in-ear sound generation requires less mechanical displacement, MEMS loudspeakers have shown interesting results in terms of frequency response and sound pressure level when evaluated in couplers [8–10]. A few MEMS loudspeakers for in-ear applications have demonstrated to fulfill the industry required performances, with electrostatic [8] and piezoelectric transduction [9,10].

The second most important characteristic of a loudspeaker after its frequency response is its linearity. This characteristic is often evaluated using the Total Harmonic Distortion (THD). The THD is the ratio of the power of all harmonics over the total power. Studies have shown that depending on the order of the harmonics, the human ear can be sensitive

<footnote n="1">¹ University Grenoble Alpes, CEA, Leti, F-38000 Grenoble, France; thierry.hilt@cea.fr (T.H.); fabrice.casset@cea.fr (F.C.); christophe.poulain@cea.fr (C.P.); gwenael.le-rhun@cea.fr (G.L.R.); franklin.pavageau@cea.fr (F.P.); hugo.kuentz@cea.fr (H.K.); mikael.colin@cea.fr (M.C.)</footnote>

<footnote n="2">² Laboratoire d'Acoustique de l'Université du Mans, LAUM-UMR 6613 CNRS, Le Mans Université, F-72085 Le Mans, France; stephane.durand@univ-lemans.fr</footnote>

<footnote>* Correspondence: Romain.Liechti@cea.fr</footnote>
to a THD as low as 0.1% [11–14]. Electrostatic transduction is inherently nonlinear, and the reachable displacement is not sufficient to advantageously replace non-MEMS electrodynamic loudspeakers. The most promising transduction is the piezoelectric transduction and more specifically the one using PZT (lead zirconate titanate) actuators, providing a high transduction factor for low actuation voltages in the case of thin films. Unfortunately, thin-film PZT shows a ferroelectric and electrostrictive behavior, creating a wide variety of nonlinearities [15–18].

To model the nonlinearities of classical electrodynamic loudspeakers, state space models, port-Hamiltonian systems, Hammerstein models and power series have been extensively used [19–25]. State space models are the most commonly used models to simulate nonlinear effects in electrodynamic loudspeaker and have shown accurate results. However, for piezoelectric transduction, such a model has not been reported yet. Recently, a model for the simulation of nonlinear effects in coupler for electrostatic transduction was reported [26]. Accurately, modeling the nonlinear behavior of loudspeakers from static nonlinear measurements of the transduction transfer function could be an efficient way for estimating the total harmonic distortion, without using expensive prototypes. In addition, an accurate model of the nonlinear behavior of the loudspeaker could help developing methods to actively reduce the nonlinearities, as it has been extensively performed with electrodynamic non-MEMS loudspeakers [27]. In addition, accurately simulating the nonlinearities generated by a loudspeaker could be used to simulate the subjective audio quality of the loudspeaker in the early phases of the design, using models of the perceptual components of the loudspeaker [28–31]. The aim of this study is to provide such a model to evaluate the distortion generated by a cantilever-based piezoelectric MEMS loudspeaker using static measurements and a lumped equivalent circuit in order to shorten the linearity optimization step of this type of loudspeakers in the design phase. Parameters are adapted to compare the results with the measurements of the Usound UT-P 2018, in order to validate the model.

## 2. Nonlinear State Space Model

As most loudspeakers behave as a simple mechanical oscillator for a wide range of frequencies, they can be well approximated as a first-order mechanical oscillator, two transformers and an electrical circuit. In the case of piezoelectric loudspeakers, the usual gyrator used for the electromechanical transduction is replaced by a transfomer, and the electrical circuit is a simple capacitor. The equivalent electrical circuit of the loudspeaker detailed in [32] is presented in Figure 1, where $R_g$ is the output resistance of the amplifier and the resistance of the connection lines, $C_p$ is the piezoelectric capacitance, $\gamma$ is the transduction ratio of the piezoelectric layer, $R_{ms}$ represents the viscous losses of the mechanical oscillator, $M_{ms}$ is the moving mass of the loudspeaker, $C_{ms}$ is the compliance of the loudspeaker and $S_d$ is the effective radiating area of the membrane.

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Equivalent electrical circuit of the loudspeaker.</figcaption></figure>

In the case of cantilever-based piezoelectric loudspeakers, the transduction factor $\gamma$ can be estimated using the unimorph piezoelectric cantilever model described in [33]. A schematic representation of a unimorph cantilever is depicted in Figure 2 with the relevant dimensions and axes. Considering the cantilever clamped on one end and guided on the
other end and using a null displacement at the tip of the cantilever, the blocked force $F_{bl}$ generated by the actuators as a function of the input voltage can be written:

$$F_{bl} = \frac{6Wt_p}{4s_{11}^p L} \frac{AB(B+1)}{AB+1} \left(1 - \frac{k_p}{k_p + k_a}\right) e_{31}(s_{11}^p + s_{12}^p)u_c \quad (1)$$

where W is the width of the actuator, $t_p$ is the thickness of the piezoelectric layer, $s_{11}^p$ and $s_{12}^p$ are the elastic compliance of the piezoelectric layer at constant electric field in the 1 and 2 direction, L is the length of the cantilever, $t_m$ is the thickness of the elastic layer, $A = s_{11}^p/s_{11}^m$, $B = t_m/t_p$, and $k_a$ and $k_p$ are the apparent stiffnesses of the active and passive parts of length L and $L_p$ of the actuators. The transduction factor $\gamma$ is the ratio of the blocked force $F_{bl}$ over the voltage $u_c$.

<figure><img src="image_7.png" /><figcaption>Figure 2. Schematic representation of a unimorph piezoelectric cantilever.</figcaption></figure>

In our case, the load impedance $Z_f$ is the acoustical impedance seen by the front side of the membrane and is equal to the acoustical impedance of the coupler. Due to the size of diaphragm of MEMS loudspeakers, which are mostly smaller than 1 cm², the back acoustic radiation impedance of the loudspeaker is not considered, due to its insignificant effect on the frequency response. A coupler is a mechanical device used to simulate the acoustical impedance of the average human ear, following the definition given in the IEC standard IEC 60318-4 [34]. The known electrical equivalent circuit is depicted in Figure 3.

<figure><img src="image_8.png" /><figcaption>Figure 3. Lumped equivalent network of the coupler.</figcaption></figure>

Using the circuit resulting of the combination of Figures 1 and 3, time domain differential equations can be derived. For the system to be fully defined, one state variable is needed for each capacitive or inductive element. The displacement $x$ is also chosen, because it is used to evaluate the apparent stiffness of the actuators. The thirteen equations of the state variables are given in Appendix A. Equations (A1)–(A13) can be used in a matrix form and solved by using a differential approximation scheme [35]. The state space equation then writes:

$$\frac{\partial x}{\partial t} = Ax + Bu_{in}(t) \quad (2)$$
with

$$ \mathbf{x} = \begin{bmatrix} u_c(t) \\ v(t) \\ x(t) \\ F_{Cms}(t) \\ p_{C1}(t) \\ Q_2(t) \\ p_{C2}(t) \\ Q_{L3}(t) \\ p_{C3}(t) \\ Q_4(t) \\ p_{C4}(t) \\ Q_{out}(t) \\ p_{out}(t) \end{bmatrix}, \quad (3) $$

$$ A = \begin{bmatrix} -\frac{1}{R_g C_p(u_c)} & -\frac{\gamma(u_c)}{C_p(u_c)} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{\gamma(u_c)}{M_{ms} + S_d^2 L_1} & -\frac{R_{ms}}{M_{ms} + S_d^2 L_1} & 0 & -\frac{1}{M_{ms} + S_d^2 L_1} & -\frac{S_d}{M_{ms} + S_d^2 L_1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{C_{ms}(x)} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{S_d}{C_1} & 0 & 0 & 0 & -\frac{1}{C_1} & 0 & -\frac{1}{C_1} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{L_2} & -\frac{R_2}{L_2} & -\frac{1}{L_2} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{C_2} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{L_3} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{C_3} & 0 & -\frac{1}{C_3} & 0 & -\frac{1}{C_3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{L_4} & -\frac{R_4}{L_4} & -\frac{1}{L_4} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{C_4} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{L_5} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{C_5} & 0 & 0 \end{bmatrix} \quad (4) $$

and

$$ B = \begin{bmatrix} \frac{1}{R_g C_p(u_c)} \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}. \quad (5) $$

Equation (2) is of the form:

$$ \dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u} \quad (6) $$

and thus, matrices can be discretized using a bilinear transform:

$$ \mathbf{A}_d = \left( \mathbf{I} - \frac{1}{2} \mathbf{A} T_s \right)^{-1} \left( \mathbf{I} + \frac{1}{2} \mathbf{A} T_s \right) \quad (7) $$
and

$$ \mathbf{B}_d = \left( \mathbf{I} - \frac{1}{2} \mathbf{A} \mathbf{T}_s \right)^{-1} \mathbf{B} \mathbf{T}_s \qquad (8) $$

with $T_s = 1/f_s$ the sampling period. Using matrices $A_d$ and $B_d$ and an input signal, the output can be computed:

$$ \mathbf{x}(k+1) = \mathbf{A}_d \mathbf{x}(k) + \mathbf{B}_d \mathbf{u}(k) \qquad (9) $$

From the computed output pressure in the state vector $\mathbf{x}$, the THD, which is the ratio between the sum of effective values of the harmonics and the effective value of the fundamental harmonic and other harmonics, is calculated using:

$$ \text{THD} = 100 \cdot \frac{\sqrt{\sum_{h=2}^{H} v_h^2}}{\sqrt{\sum_{h=1}^{H} v_h^2}} \qquad (10) $$

### 3. Static Measurements

As other loudspeakers state space models, this model uses nonlinear functions as input parameters. The nonlinear functions can be measured dynamically using adaptive state space models, calculated analytically, computed with finite element models or determined from static measurement. In this paper, we propose to use static measurements of the nonlinear transfer functions, measured on a commercial loudspeaker [36] and on PZT test devices. Since internal dimensions of the loudspeaker are unknown and difficult to measure, the lumped element parameters of the equivalent circuit $C_{ms}$, $S_d$, $M_{ms}$ and $C_p$ are taken directly from the data sheet.

#### 3.1. Nonlinear Electromechanical Function

The $e_{31,f}$ piezoelectric coefficient of thin-film PZT is extracted from the $d_{33,f}$ measurement, with the procedure described in [37], as a function of the applied electric field, using the DBLI, on a standard PZT test device. As the piezoelectric coefficient $e_{31}$ is calculated for a piezoelectric material on wafer with this specific machine, it would be difficult to measure it on the loudspeaker directly. The parameter $e_{31}$ is known to be highly dependent on the applied electric field, due to intrinsic and extrinsic various effects, such as electrostriction, ferroelectric behavior and the domain movements. With the complete physical interactions being not completely understood, these nonlinearities are often simulated using phenomenological models, whose parameters are deduced from measurements [38,39]. The parameter $e_{31,f}$ and the fitted logarithmic function are depicted in Figure 4.

<figure><img src="image_9.png" /><figcaption>Figure 4. Measured and fitted function for the nonlinear transverse piezoelectric coefficient.</figcaption></figure>
As it has been reported in multiple papers, the transversal piezoelectric coefficient of piezoelectric thin film increases for low electric fields, until it converges to the maximum value, due to the alignment of the domains of the material [40–44]. After the saturation of the rotation of the domains, the piezoelectric coefficient decreases slowly. The fitted function of $e_{31,f}$ is then used in Equation (1) to calculate the nonlinear transduction factor $\gamma(u_c)$. Depending on the actuators geometry, the fitted $e_{31,f}$ function can be used in any other electromechanical transduction function.

## 3.2. Nonlinear Stiffness

Using the setup depicted in Figure 5, the nonlinear relation between the force and the displacement of the commercial loudspeaker, not subjected to any voltage, is measured.

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Schematic representation of a nanoindenter.</figcaption></figure>

The commercial loudspeaker is placed on the sample loader. The nanoindenter, equipped with a spherical tip of radius 50 µm, applies a precise and well-known force on the moving part of the sample, until reaching the previously defined maximum displacement of 30 µm. After a processing of the raw data by the embedded software, the curve of the force applied on the sample and of the displacement of the shaft is retrieved and presented in Figure 6.

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Nonlinear stiffness of the loudspeaker.</figcaption></figure>

A third-order polynomial function, known to accurately represent the nonlinear behavior of clamped-guided structures, is used as the nonlinear fitting function. The linear
approximation is calculated from the $V_{AS}$, which is the equivalent compliance volume value given in the loudspeaker data sheet, using:

$$C_{ms} = \frac{V_{AS}}{\rho_0 c_0^2 S_d^2} \qquad (11)$$

with $\rho_0$ the density of air at room temperature, $c_0$ the speed of sound in air at room temperature and $S_d$ the effective radiating surface of the loudspeaker [45]. From the fitted function, the nonlinear compliance $C_{ms}(x)$ can be calculated. This nonlinear function can also be calculated knowing the geometry of the loudspeaker and a finite element model or analytically.

### 3.3. Capacitance

As the PZT thin film is ferroelectric, the capacitance created by the piezoelectric layer between the electrodes varies with the applied electric field. The nonlinear relation was measured on a layer of PZT using the AixACCT TF Analyzer 2000. The measured variation of the capacitance and the fitted Lorentz function are presented in Figure 7.

<figure><img src="image_12.png" /><figcaption>Figure 7. Non linear capacitance of the loudspeaker.</figcaption></figure>

From the fitted Lorentzian function, the nonlinear capacitance $C_p(u_c)$ can be calculated. The curves can be adjusted in order to match the value of the data sheet of the loudspeaker.

## 4. Results

Using Equation (9) and a scaled dirac whose amplitude at the first sample is $f_s u_p/2$ as an input signal, the linear frequency response of the loudspeaker inside the coupler is computed. The computed frequency response is compared to the frequency response of the commercial loudspeaker, in Figure 8. The parameters for the circuits presented on Figures 1 and 3 are derived or directly taken from the data sheet of the loudspeaker [36] and taken from the literature for the coupler [46]. Although the influence of $R_g$ is negligible, it is set to 1 mΩ to have a properly scaled $A_d$ matrix. The value of $M_{ms}$ is computed from the previously determined value of $C_{ms}$ and the value of the resonance frequency given in the data sheet using:

$$M_{ms} = \frac{1}{4\pi^2 f_s^2 C_{ms}} \qquad (12)$$

The value of $R_{ms}$ is determined using the quality factor of the data sheet using:

$$R_{ms} = \frac{1}{Q} \sqrt{\frac{M_{ms}}{C_{ms}}} \qquad (13)$$
All parameters are summarized in Table 1.

<table><caption>**Table 1.** Table of lumped elements circuits parameters.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>R<sub>g</sub></td><td>1</td><td>mΩ</td></tr><tr><td>C<sub>p</sub></td><td>39</td><td>nF</td></tr><tr><td>R<sub>ms</sub></td><td>5 × 10<sup>-3</sup></td><td>Ns/m</td></tr><tr><td>M<sub>ms</sub></td><td>1.16</td><td>mg</td></tr><tr><td>C<sub>ms</sub></td><td>3 × 10<sup>-3</sup></td><td>m/N</td></tr><tr><td>S<sub>d</sub></td><td>12</td><td>mm<sup>2</sup></td></tr><tr><td>L<sub>1</sub></td><td>82.9</td><td>Pas<sup>2</sup>/m<sup>3</sup></td></tr><tr><td>L<sub>2</sub></td><td>9400</td><td>Pas<sup>2</sup>/m<sup>3</sup></td></tr><tr><td>L<sub>3</sub></td><td>130.3</td><td>Pas<sup>2</sup>/m<sup>3</sup></td></tr><tr><td>L<sub>4</sub></td><td>983.8</td><td>Pas<sup>2</sup>/m<sup>3</sup></td></tr><tr><td>L<sub>5</sub></td><td>133.4</td><td>Pas<sup>2</sup>/m<sup>3</sup></td></tr><tr><td>C<sub>1</sub></td><td>1 × 10<sup>-12</sup></td><td>Pa/m<sup>3</sup></td></tr><tr><td>C<sub>2</sub></td><td>1.9 × 10<sup>-12</sup></td><td>Pa/m<sup>3</sup></td></tr><tr><td>C<sub>3</sub></td><td>1.5 × 10<sup>-12</sup></td><td>Pa/m<sup>3</sup></td></tr><tr><td>C<sub>4</sub></td><td>2.1 × 10<sup>-12</sup></td><td>Pa/m<sup>3</sup></td></tr><tr><td>C<sub>5</sub></td><td>1.517 × 10<sup>-12</sup></td><td>Pa/m<sup>3</sup></td></tr><tr><td>R<sub>1</sub></td><td>50.6 × 10<sup>6</sup></td><td>Pas/m<sup>3</sup></td></tr><tr><td>R<sub>2</sub></td><td>31.1 × 10<sup>6</sup></td><td>Pas/m<sup>3</sup></td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 8.** Simulated frequency response of the commercial loudspeaker and measured frequency response of the loudspeaker adapted from [36].</figcaption></figure>

Despite some high-frequency mismatches, due to the modal behavior of the membrane above 20 kHz, and a slight difference between the quality factors of the coupler resonances at, respectively, 12,710 and 20,970 Hz, the simulated frequency response gives a good approximation of the manufacturer measurement.

Using the identified nonlinear functions for parameters $\gamma$, $K_{ms}$ and $C_p$ and a pure sine of 1 kHz of amplitude 1.41 V as an input signal, the Fourier transform of the output signal can be calculated and is presented in Figure 9, for different input signals, as well as the linear frequency response.
<figure><img src="image_14.png" /><figcaption>(a) 100 Hz</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) 1 kHz</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) Time domain 1 kHz</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) 70 Hz + 600 Hz</figcaption></figure>

**Figure 9.** Nonlinear loudspeaker model responses to different input signals with the linear frequency response as a reference in frequency and time domain.

In Figures 9a,b, the response of the loudspeaker to a pure sine of, respectively, 100 Hz and 1 kHz are given, and one can see the second and third harmonic clearly appear. In Figure 9c, the displacement and input voltage in time domain are given. In Figure 9d, the response of the loudspeaker to a sum of two pure sines with $f_1 = 70$ Hz and $f_2 = 600$ Hz, is depicted, and one can see that the intermodulation appears with the components of frequencies $f_2 \pm f_1$, $f_2 \pm 2f_1$, $2f_2 \pm f_1$ and $3f_2 \pm f_1$.

The THD is computed point by point using a pure sine as the input signal for 200 points. The simulated THD is compared to the THD measured by the manufacturer in Figure 10. The pressure signal is low pass filtered at 16 kHz before the THD estimation, as it is described in the data sheet of the loudspeaker. There is a good agreement between the absolute level of the THD between 40 and 200 Hz, where the difference between the two curves is below 20%. Above 500 Hz, the difference between the simulated and measured frequency responses is the cause of the difference between the THD curves and around the distortion peak at about 6.4 kHz, the difference is maximum and above 100%. Qualitatively, the shape of the curve around the resonance of the loudspeaker has a similar shape, capturing the effect of the resonance of the loudspeaker and the one of the coupler on the THD.
<figure><img src="image_18.png" /><figcaption>**Figure 10.** Simulated THD of the commercial loudspeaker and measured THD of the loudspeaker adapted from [36].</figcaption></figure>

In Figure 11a, the respective contributions of $\gamma(u_c)$, $C_{ms}(x)$ and $C_p(u_c)$ to the total harmonic distortion are shown. As the displacement, the harmonic distortion generated by $C_{ms}(x)$ is maximum below the resonance frequency of the loudspeaker and negligible above the resonance frequency. The harmonic distortion generated by the transduction factor $\gamma(u_c)$ is modulated by the frequency response of the loudspeaker but remains high in the whole audio bandwidth. The contribution of $C_p(u_c)$ is three orders of magnitude below the other ones and then negligible. In Figure 11b, the THD for three different input levels is presented.

<figure><img src="image_19.png" /><figcaption>**Figure 11.** THD generated by individual parameters (a) and THD as a function of input level (b).</figcaption></figure>

**5. Conclusions**

The model presented in this paper allows one to evaluate the THD generated by a cantilever-based MEMS loudspeaker, by only measuring the nonlinear relations between the stiffness and the displacement, the capacitance and the electric field, and the force and the electric field. It also could allow one to use predistortion techniques to reduce the said total harmonic distortion in MEMS loudspeakers. Although open-loop correction of loudspeaker nonlinearities is tedious for non-MEMS loudspeakers, due to large manufacturing discrepancies, this technique could be more suited for MEMS design, due to the tightest tolerances of the manufacturing processes. This model is also very efficient in terms of computational power compared to a finite element modeling solution and could be used in the early stages of piezoelectric MEMS loudspeakers development, to evaluate
the linear performance of the loudspeaker. The discrepancies between the model and the measurement can be attributed to the assumptions made in the model of the loudspeaker. For example, the nonlinear behavior of the PZT is known to widely vary with frequency and to create hysteresis. In addition, other different nonlinear parameters were not taken into account, such as the nonlinearities due to high velocity in the coupler slits, but can be added easily to the model.

**Author Contributions:** Conceptualization, R.L., T.H., F.C. and S.D.; methodology, R.L.; validation, R.L.; formal analysis, R.L.; investigation, C.P., F.P., H.K., G.L.R. and R.L.; data curation, R.L.; writing—original draft preparation, R.L.; writing—review and editing, R.L., T.H., F.C., S.D., G.L.R. and C.P.; supervision, T.H., F.C., S.D. and M.C.; project administration, M.C.; funding acquisition, T.H. and M.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable

**Informed Consent Statement:** Not applicable

**Conflicts of Interest:** The authors declare no conflict of interest.

## Abbreviations

The following abbreviations are used in this manuscript:

THD Total Harmonic Distortion

MEMS Micro Electromechanical Systems

IEC International Electrotechnical Commission

DUT Device Under Test

DBLI Double-Beam Laser Interferometer

PZT Lead Zirconate Titanate

## Appendix A. State Space Equations

$$ \frac{\partial u_c(t)}{\partial t} = \frac{1}{R_g C_p(u_c)} u_{in}(t) - \frac{1}{R_g C_p(u_c)} u_c(t) - \frac{\gamma(u_c)}{C_p(u_c)} v(t) \quad (A1) $$

$$ \frac{\partial v(t)}{\partial t} = \frac{\gamma(u_c)}{M_{ms} + S_d^2 L_1} u_c(t) - \frac{R_{ms}}{M_{ms} + S_d^2 L_1} v(t) - \frac{1}{M_{ms} + S_d^2 L_1} F_{C_{ms}}(t) - \frac{S_d}{M_{ms} + S_d^2 L_1} p_{C_1}(t) \quad (A2) $$

$$ \frac{\partial F_{C_{ms}}}{\partial t} = \frac{1}{C_{ms}(x)} v(t), \quad (A3) $$

$$ \frac{\partial x(t)}{\partial t} = v(t), \quad (A4) $$

$$ \frac{\partial p_{C_1}(t)}{\partial t} = \frac{S_d}{C_1} v(t) - \frac{1}{C_1} Q_2(t) - \frac{1}{C_1} Q_{L_3}(t) \quad (A5) $$

$$ \frac{\partial Q_2(t)}{\partial t} = \frac{1}{L_2} p_{C_1}(t) - \frac{1}{L_2} p_{C_2}(t) - \frac{R_2}{L_2} Q_2(t) \quad (A6) $$

$$ \frac{\partial p_{C_2}(t)}{\partial t} = \frac{1}{C_2} Q_2(t) \quad (A7) $$

$$ \frac{\partial Q_{L_3}(t)}{\partial t} = \frac{1}{L_3} p_{C_1}(t) - \frac{1}{L_3} p_{C_3}(t) \quad (A8) $$
$$ \frac{\partial p_{C_3}(t)}{\partial t} = \frac{1}{C_3} Q_{L_3}(t) - \frac{1}{C_3} Q_4(t) - \frac{1}{C_3} Q_{out}(t) \qquad (A9) $$

$$ \frac{\partial Q_4(t)}{\partial t} = \frac{1}{L_4} p_{C_3}(t) - \frac{1}{L_4} Q_4(t) - \frac{1}{L_4} p_{C_4}(t) \qquad (A10) $$

$$ \frac{\partial p_{C_4}(t)}{\partial t} = \frac{1}{C_4} Q_4(t) \qquad (A11) $$

$$ \frac{\partial Q_{out}(t)}{\partial t} = \frac{1}{L_5} p_{C_3}(t) - \frac{1}{L_5} p_{out}(t) \qquad (A12) $$

$$ \frac{\partial p_{out}(t)}{\partial t} = \frac{1}{C_5} Q_{out}(t) \qquad (A13) $$

## References

1. Shahosseini, I.; Lefeuvre, E.; Moulin, J.; Martincic, E.; Woytasik, M.; Lemarquand, G. Optimization and Microfabrication of High Performance Silicon-Based MEMS Microspeaker. *IEEE Sens. J.* **2013**, *13*, 273–284. [CrossRef]

2. Shahosseini, I. Vers des Micro-Haut-Parleurs à Hautes Performances Électroacoustiques en Technologie Silicium. Ph.D. Thesis, Université Paris Sud, Paris, France, 2012; p. 206.

3. Wang, H.; Feng, P.X.L.; Xie, H. A Dual-Electrode MEMS Speaker Based on Ceramic PZT with Improved Sound Pressure Level by Phase Tuning. In Proceedings of the 2021 IEEE 34th International Conference on Micro Electro Mechanical Systems (MEMS), Gainesville, FL, USA, 25–29 January 2021; pp. 701–704. [CrossRef]

4. Wang, Q.; Yi, Z.; Ruan, T.; Xu, Q.; Yang, B.; Liu, J. Obtaining High SPL Piezoelectric MEMS Speaker via a Rigid-Flexible Vibration Coupling Mechanism. *J. Microelectromech. Syst.* **2021**, *30*, 725–732. [CrossRef]

5. Ko, S.C.; Kim, Y.C.; Lee, S.S.; Choi, S.H.; Kim, S.R. Micromachined piezoelectric membrane acoustic device. *Sens. Actuators A Phys.* **2003**, *103*, 130–134. [CrossRef]

6. Wang, H.; Chen, Z.; Xie, H. A high-SPL piezoelectric MEMS loud speaker based on thin ceramic PZT. *Sens. Actuators A Phys.* **2020**, *309*, 112018. [CrossRef]

7. Wang, H.; Li, M.; Yu, Y.; Chen, Z.; Ding, Y.; Jiang, H.; Xie, H. A Piezoelectric MEMS Loud Speaker Based on Ceramic PZT. In Proceedings of the 2019 20th International Conference on Solid-State Sensors, Actuators and Microsystems & Eurosensors XXXIII (Transducers & Eurosensors XXXIII), Berlin, Germany, 23–27 June 2019; pp. 857–860. [CrossRef]

8. Kaiser, B.; Langa, S.; Ehrig, L.; Stolz, M.; Schenk, H.; Conrad, H.; Schenk, H.; Schimmanz, K.; Schuffenhauer, D. Concept and proof for an all-silicon MEMS micro speaker utilizing air chambers. *Microsyst. Nanoeng.* **2019**, *5*, 43. [CrossRef]

9. Stoppel, F.; Eisermann, C.; Gu-Stoppel, S.; Kaden, D.; Giese, T.; Wagner, B. Novel membrane-less two-way MEMS loudspeaker based on piezoelectric dual-concentric actuators. In Proceedings of the 2017 19th International Conference on Solid-State Sensors, Actuators and Microsystems (Transducers), Kaohsiung, Taiwan, 18–22 June 2017; pp. 2047–2050. [CrossRef]

10. Stoppel, F.; Mannchen, A.; Niekiel, F.; Beer, D.; Giese, T.; Wagner, B. New integrated full-range MEMS speaker for in-ear applications. In Proceedings of the 2018 IEEE Micro Electro Mechanical Systems (MEMS), Belfast, UK, 21–25 January 2018; pp. 1068–1071. [CrossRef]

11. Fielder, L.; Benjamin, E. Subwoofer Performance for Accurate Reproduction of Music. *J. Audio Eng. Soc.* **1988**, *36*, 443–456.

12. Vanderkooy, J.; Krauel K.B. Another View of Distortion Perception. In *Audio Engineering Society Convention 133*; Audio Engineering Society: New York, NY, USA, 2012; p. 4.

13. Larson, J.; DellaSala, G. A Tutorial on the Audibility of Loudspeaker Distortion at Bass Frequencies. In *Audio Engineering Society Convention 143*; Audio Engineering Society: New York, NY, USA, 2017; p. 10.

14. Temme, S.; Brunet, P.; Qarabaqi, P. Measurement of Harmonic Distortion Audibility Using a Simplified Psychoacoustic Model-Updated. In *Audio Engineering Society Conference: 51st International Conference: Loudspeakers and Headphones*; Audio Engineering Society: New York, NY, USA, 2013; p. 9.

15. Damjanovic, D. Hysteresis in Piezoelectric and Ferroelectric Materials. In *The Science of Hysteresis*; Elsevier: Amsterdam, The Netherlands, 2006; pp. 337–465. [CrossRef]

16. Jayendran, R.; Arockiarajan, A. Non-linear electromechanical response of 1–3 type piezocomposites. *Int. J. Solids Struct.* **2013**, *50*, 2259–2270. [CrossRef]

17. Royston, T.J.; Houston, B.H. Modeling and measurement of nonlinear dynamic behavior in piezoelectric ceramics with application to 1–3 composites. *J. Acoust. Soc. Am.* **1998**, *104*, 2814–2827. [CrossRef]
18. Lucke, P.; Bayraktar, M.; Birkhölzer, Y.A.; Nematollahi, M.; Yakshin, A.; Rijnders, G.; Bijkerk, F.; Houwman, E.P. Hysteresis, Loss and Nonlinearity in Epitaxial PbZr<sub>0.55</sub>Ti<sub>0.45</sub>O<sub>3</sub> Films: A Polarization Rotation Model. *Adv. Funct. Mater.* **2020**, *30*, 2005397. [CrossRef]

19. Rébillat, M.; Hennequin, R.; Corteel, E.; Katz, B.F. Identification of cascade of Hammerstein models for the description of nonlinearities in vibrating devices. *J. Sound Vib.* **2011**, *330*, 1018–1038. [CrossRef]

20. Falaize, A.; Hélie, T. Passive modelling of the electrodynamic loudspeaker: From the Thiele–Small model to nonlinear port-Hamiltonian systems. *Acta Acust. 2020*, *4*, 1. [CrossRef]

21. Falaize, A.; Papazoglou, N.; Hélie, T.; Lopes, N. Compensation of loudspeaker's nonlinearities based on flatness and port-Hamiltonian approach. In *22ème Congrès Français de Mécanique*, Association Française de Mécanique: Lyon, France, 2015; p. 11.

22. Novak, A.; Simon, L.; Lotton, P.; Merit, B.; Gilbert, J. Nonlinear Analysis and Modeling of Electrodynamic Loudspeakers. In *10ème Congrès Français d'Acoustique*, Société Française d'Acoustique: Paris, France, 2010; p. 6.

23. Klippel, W. Large Signal Performance of Tweeters, Micro Speakers and Horn Drivers. In *Audio Engineering Society Convention 118*; Audio Engineering Society: New York, NY, USA, 2005; p. 19.

24. Klippel, W. Loudspeaker Nonlinearities–Causes, Parameters, Symptoms. In *Audio Engineering Society Convention 119*; Audio Engineering Society: New York, NY, USA, 2005; p. 36.

25. King, A.; Agerkvist, F. Fractional Derivative Loudspeaker Models for Nonlinear Suspensions and Voice Coils. *J. Audio Eng. Soc.* **2018**, *66*, 525–536. [CrossRef]

26. Monsalve, J.M.; Melnikov, A.; Kaiser, B.; Schuffenhauer, D.; Stolz, M.; Ehrig, L.; Schenk, H.A.G.; Conrad, H.; Schenk, H. Large-Signal Equivalent-Circuit Model of Asymmetric Electrostatic Transducers. *IEEE/ASME Trans. Mechatron.* **2021**, 1–11. [CrossRef]

27. Klippel, W. Optimal Design of Loudspeakers with Non-Linear Control. In *Audio Engineering Society Conference: 32nd International Conference: DSP For Loudspeakers*; Audio Engineering Society: New York, NY, USA, 2007; p. 11.

28. Moore, B.C.J.; Tan, C.T. Perceived naturalness of spectrally distorted speech and music. *J. Acoust. Soc. Am.* **2003**, *114*, 408–419. [CrossRef]

29. Tan, C.T.; Moore, B.C.J. The Effect of Nonlinear Distortion on the Perceived Quality of Music and Speech Signals. *J. Audio Eng. Soc.* **2003**, *51*, 20.

30. Moore, B.C.J. Computational models for predicting sound quality. *Acoust. Sci. Technol.* **2020**, *41*, 75–82. [CrossRef]

31. Olsen, S.L.; Agerkvist, F.; MacDonald, E.; Stegenborg-Andersen, T.; Volk, C.P. Modelling the Perceptual Components of Loudspeaker Distortion. In *Audio Engineering Society Convention 140*; Audio Engineering Society: New York, NY, USA, 2016; p. 9.

32. Liechti, R.; Durand, S.; Hilt, T.; Casset, F.; Dieppedale, C.; Verdot, T.; Colin, M. A Piezoelectric MEMS Loudspeaker Lumped and FEM models. In Proceedings of the 2021 22nd International Conference on Thermal, Mechanical and Multi-Physics Simulation and Experiments in Microelectronics and Microsystems (EuroSimE), Online, 19–21 April 2021; pp. 1–8. [CrossRef]

33. Wang, Q.M.; Du, X.H.; Xu, B.; Cross, L.E. Electromechanical coupling and output efficiency of piezoelectric bending actuators. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **1999**, *46*, 638–646. [CrossRef]

34. International Electrotechnical Commission. *Electroacoustics: Simulators of Human Head and Ear. Part 4, Partie 4*; OCLC: 954221482; International Electrotechnical Commission: Geneva, Switzerland, 2010.

35. Fadali, M.S.; Visioli, A. *Digital Control Engineering: Analysis and Design*, 2nd ed.; Academic Press, Elsevier: Amsterdam, The Netherlands, 2013.

36. Usound. Achelous UT-P 2018 Datasheet. 2020. Available online: https://www.usound.com/wp-content/uploads/2020/01/2018_Achelous-UT-P-2018-Datasheet.pdf (accessed on 23 November 2021).

37. Sivaramakrishnan, S.; Mardilovich, P.; Schmitz-Kempen, T.; Tiedke, S. Concurrent wafer-level measurement of longitudinal and transverse effective piezoelectric coefficients ($d_{33,f}$ and $e_{31,f}$) by double beam laser interferometry. *J. Appl. Phys.* **2018**, *123*, 014103. [CrossRef]

38. Marton, P.; Rychetsky, I.; Hlinka, J. Domain walls of ferroelectric BaTiO₃ within the Ginzburg-Landau-Devonshire phenomenological model. *Phys. Rev. B* **2010**, *81*, 144125. [CrossRef]

39. Mehling, V.; Tsakmakis, C.; Gross, D. Phenomenological model for the macroscopical material behavior of ferroelectric ceramics. *J. Mech. Phys. Solids* **2007**, *55*, 2106–2141. [CrossRef]

40. Tsujiura, Y.; Kawabe, S.; Kurokawa, F.; Hida, H.; Kanno, I. Comparison of effective transverse piezoelectric coefficients $e_{31,f}$ of Pb(Zr,Ti)O₃ thin films between direct and converse piezoelectric effects. *Jpn. J. Appl. Phys.* **2015**, *54*, 10NA04. [CrossRef]

41. Völker, B.; Kamlah, M.; Wang, J. Phase-field Modeling of Ferroelectric Materials. In Proceedings of the COMSOL Conference 2009, Milanop, Italy, 14–16 October 2009; p. 28.

42. Tan, G.; Maruyama, K.; Kanamitsu, Y.; Nishioka, S.; Ozaki, T.; Umegaki, T.; Hida, H.; Kanno, I. Crystallographic contributions to piezoelectric properties in PZT thin films. *Sci. Rep.* **2019**, *9*, 7309. [CrossRef] [PubMed]

43. Wang, Y.; Cheng, H.; Yan, J.; Chen, N.; Yan, P.; Ouyang, J. Nonlinear electric field dependence of the transverse piezoelectric response in a (001) ferroelectric film. *Scr. Mater.* **2020**, *189*, 84–88. [CrossRef]

44. Le Rhun, G.; Dieppedale, C.; Wague, B.; Querne, C.; Enyedi, G.; Perreau, P.; Montmeat, P.; Licitra, C.; Fanget, S. Transparent PZT MIM Capacitors on Glass for Piezoelectric Transducer Applications. In Proceedings of the 2019 20th International Conference on Solid-State Sensors, Actuators and Microsystems & Eurosensors XXXIII (Transducers & Eurosensors XXXIII), Berlin, Germany, 23–27 June 2019; pp. 1800–1802. [CrossRef]
45. Kleiner, M. *Electroacoustics*; CRC Press: Boca Raton, FL, USA, 2013; p. 620.

46. Jönsson, S.; Liu, B.; Nielsen, L.B.; Schuhmacher, A. Simulation of Couplers. In *Audio Engineering Society Workshop* 7; Audio Engineering Society: New York, NY, USA, 2003; p. 29.