import matplotlib.pyplot as plt
Funtional Transformation Method Utilities
This notebook contains a set of utility functions for the functional transformation method.
::: {#cell-4 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’ 6=‘i’}
from dataclasses import asdict, dataclass
import numpy as np
from opt_einsum import contract as einsum
from scipy.integrate import simpson, trapezoid
from scipy.signal import cont2discrete
from scipy.special import jn_zeros, jv
from tabulate import tabulate
:::
::: {#cell-5 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
@dataclass
class PhysicalParameters:
def asdict(self):
return asdict(self)
def tabulate(self) -> str:
dict[str, str] = {
param_symbols: "A": "$A$",
"I": "$I$",
"rho": "$\\rho$",
"E": "$E$",
"d1": "$d_1$",
"d3": "$d_3$",
"Ts0": "$T_{s0}$",
"length": "$\\ell$",
"h": "$h$",
"l1": "$l_1$",
"l2": "$l_2$",
"nu": "$\\nu$",
"d0": "$d_0$",
"d2": "$d_2$",
"Tm": "$T_m$",
"f_1": "$f_1$",
"r0": "$r_0$",
"f0": "$f_0$",
}
dict[str, str] = {
param_units: "A": "$m^2$",
"I": "$m^4$",
"rho": "$kg/m^3$",
"E": "$Pa$",
"d1": "$kg/(ms)$",
"d3": "$kg\\cdot m/s$",
"Ts0": "$N$",
"length": "$m$",
"h": "$m$",
"l1": "$m$",
"l2": "$m$",
"nu": "dimensionless",
"d0": "$kg/(m^2 s)$",
"d2": "$kg/s$",
"Tm": "$N/m$",
"f_1": "$Hz$",
"r0": "$m$",
"f0": "$Hz$",
}
= []
table_data
for field_name, field_value in self.asdict().items():
if field_name in param_symbols:
= param_symbols[field_name]
symbol = param_units.get(field_name, "")
unit
table_data.append([symbol, field_value, unit])
return tabulate(
table_data,=["Parameter", "Value", "Units"],
headers="github",
tablefmt
)
@dataclass
class StringParameters(PhysicalParameters):
"""
Dataclass to store the parameters of the string. Based on a nylon B string
of the guitar. From Table 4.1 of Digital Sound Synthesis using the FTM.
Moment of intertia and damping are changed however.
"""
# fmt: off
float = 0.5188e-6 # m**2 Cross section area
A: float = 0.141e-12 # m**4 Moment of intertia
I: float = 1140 # kg/m**3 Density
rho: int = 5.4e9 # Pa Young's modulus
E: float= 8e-5 # kg/(ms) Frequency independent loss
d1: float = 1.4e-5 # kg m/s Frequency dependent loss
d3: float = 60.97 # N Tension
Ts0: float = 0.65 # m Length of the string
length: # fmt: on
@staticmethod
def piano_string():
"""
From Table 5.1 of Digital Sound Synthesis using the FTM
"""
return StringParameters(
=1.54e-6,
A=4.12e-12,
I=57.0e3,
rho=19.5e9,
E=3e-3,
d1=2e-5,
d3=2104,
Ts0=1.08,
length
)
@staticmethod
def bass_string():
"""
From Table 5.1 of Digital Sound Synthesis using the FTM
"""
return StringParameters(
=2.4e-6,
A=0.916e-12,
I=6300,
rho=5e9,
E=6e-3,
d1=1e-3,
d3=114,
Ts0=1.05,
length
)
@staticmethod
def guitar_string_D():
"""
From Table 8.2 of Simulation of Distributed Parameter Systems by Transfer Function Models
"""
return StringParameters(
=7.96e-7,
A=0.171e-12,
I=1140,
rho=5.4e9,
E=8e-5,
d1=1.4e-5,
d3=13.35,
Ts0=0.65,
length
)
@staticmethod
def guitar_string_B_schafer():
"""
From A String In a Room: Mixed-Dimensional Transfer Function Models for Sound Synthesis
"""
return StringParameters(
=0.5e-6,
A=0.17e-12,
I=1140,
rho=5.4e9,
E=8e-5,
d1=1.4e-5,
d3=60.97,
Ts0=0.65,
length
)
@property
def area_density(self):
return self.rho * self.A
@dataclass
class PlateParameters(PhysicalParameters):
"""
From Digital sound synthesis of string instruments with the functional transformation method Table 5.2.
"""
# fmt: off
float = 1.2e-3 # m Thickness
h: float = 1.08 # m Width
l1: float = 1.08 # m Height
l2: float = 144e-12 # m**4 Moment of intertia
I: float = 7.8e3 # kg/m**3 Density
rho: int = 220e9 # Pa Young's modulus
E: float = 0.29 # Poisson's ratio
nu: float= 3.2e-2 # kg/(m**2 s) Frequency independent loss
d0: float = 1.3e-3 # kg/s Frequency dependent loss
d2: float = 2010 # N/m Surface Tension
Tm: float = 10.27 # Hz Fundamental frequency
f_1: # fmt: on
@property
def surface_density(self):
return self.rho * self.h
@dataclass
class CircularDrumHeadParameters(PhysicalParameters):
"""
Kettle drum head, from Digital sound synthesis of string instruments with the functional transformation method Table 5.2.
"""
# fmt: off
float = 1.9e-4 # m Thickness
h: float = 0.328 # m Radius
r0: float = 0.57e-12 # m**4 Moment of intertia
I: float = 1.38e3 # kg/m**3 Density
rho: int = 3.5e9 # Pa Young's modulus
E: float = 0.35 # Poisson's ratio
nu: float= 0.14 # kg/(m**2 s) Frequency independent loss
d0: float = 0.32 # kg/s Frequency dependent loss
d2: float = 3990 # N/m Surface Tension
Tm: float = 143.95 # Hz Fundamental frequency
f0: # fmt: on
@property
def surface_density(self):
return self.rho * self.h
@staticmethod
def avanzini():
"""
From Section VI of "A Modular Physically Based Approach to the Sound Synthesis of Membrane Percussion Instruments"
"""
return CircularDrumHeadParameters(
=2e-4, # 0.2mmm
h=0.20, # 20cm
r0=1350.0, # 0.27 / h
rho=3.5e9, # 3.5e9
E=0.2,
nu=1.25,
d0=5e-4,
d2=1500,
Tm )
:::
::: {#cell-6 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def string_eigenvalues_sqrt(
int,
n_max_modes: float,
length:
):r"""
Compute the eigenvalues of the string.
The eigenvalues of the string are given by:
$$
\lambda_{\mu} = \left(\frac{\mu \pi}{L}\right)
$$
where $\mu$ is the mode number and $L$ is the length of the string.
"""
= np.arange(1, n_max_modes + 1)
mu return np.pi * mu / length
def string_eigenfunctions(
wavenumbers: np.ndarray,
grid: np.ndarray,-> np.ndarray:
) r"""
Compute the modes of the string.
The modes of the string are given by:
$$
K = \sin(\pi x k)
$$
where $k$ is the wavenumber and $x$ is the grid positions.
Parameters
----------
wavenumbers : np.ndarray
The wavenumbers of the string.
grid : np.ndarray
The grid positions of the string where to compute the modes.
Returns
-------
np.ndarray
The modes of the string at the given grid positions.
"""
return np.sin(np.outer(wavenumbers, grid))
def othonormal_eigenfunctions(
eigenfunctions: np.ndarray,-> np.ndarray:
) r"""
Normalize the eigenfunctions of the string.
The orthonormal eigenfunctions of the string are given by:
$$
\sqrt{\frac{2}{N}} \sin(\pi x k)
$$
where $k$ is the wavenumber and $x$ is the grid positions.
Parameters
----------
eigenfunctions : np.ndarray
The eigenfunctions of the string.
Returns
-------
np.ndarray
The orthonormal eigenfunctions of the string.
"""
return np.sqrt(2 / eigenfunctions.shape[1]) * eigenfunctions
:::
::: {#cell-7 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def plate_wavenumbers(
int,
n_max_modes_x: int,
n_max_modes_y: float,
l1: float,
l2: -> tuple[np.ndarray, np.ndarray]:
) """
Compute the wavenumbers of the plate.
Parameters:
----------
n_max_modes_x: int
The number of modes in the x direction.
n_max_modes_y: int
The number of modes in the y direction.
l1: float
The width of the plate.
l2: float
The height of the plate.
Returns:
-------
wn_x: np.ndarray
The wavenumbers in the x direction.
wn_y: np.ndarray
The wavenumbers in the y direction
"""
= np.arange(1, n_max_modes_x + 1)
mu_x = np.arange(1, n_max_modes_y + 1)
mu_y = mu_x * np.pi / l1
wavenumbers_x = mu_y * np.pi / l2
wavenumbers_y return wavenumbers_x, wavenumbers_y
def plate_eigenvalues(
# (n_max_modes_x,)
wavenumbers_x: np.ndarray, # (n_max_modes_y,)
wavenumbers_y: np.ndarray, bool = False,
squared:
):r"""
Compute the eigenvalues of the plate.
The eigenvalues of the plate are given by:
$$
\lambda_{\mu, \nu} = \pi \sqrt{\left(\frac{\mu}{L_1}\right)^2 + \left(\frac{\nu}{L_2}\right)^2}
$$
where $\mu$ and $\nu$ are the mode numbers and $L_1$ and $L_2$ are the width and height of the plate.
"""
= np.meshgrid(wavenumbers_x, wavenumbers_y)
wn_x, wn_y
# this is important because we can't simply take the square root of the sum of squares
# wn_x**4 + wn_y**4 != (wn_x**2 + wn_y**2) ** 2
return wn_x**4 + wn_y**4 if squared else wn_x**2 + wn_y**2
def plate_eigenfunctions(
# (n_max_modes_x,)
wavenumbers_x: np.ndarray, # (n_max_modes_y,)
wavenumbers_y: np.ndarray, # (n_gridpoints_x,)
x: np.ndarray, # (n_gridpoints_y,)
y: np.ndarray, -> np.ndarray:
) r"""
Compute the modes of the plate.
The modes of the plate are given by:
$$
K = \sin(\pi x k) \sin(\pi y k)
$$
where $k$ is the wavenumber and $x$ and $y$ are the grid positions.
"""
# Compute the sine values using broadcasting
= np.sin(wavenumbers_x[:, None] * x[None, :])
sin_wx_grid_x = np.sin(wavenumbers_y[:, None] * y[None, :])
sin_wy_grid_y
# Use einsum to compute the outer product and obtain the modes
return einsum("m x, n y -> m n x y", sin_wx_grid_x, sin_wy_grid_y)
:::
::: {#cell-9 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def drumhead_wavenumbers(
int,
n_max_modes: int,
m_max_modes: float,
radius: -> np.ndarray:
) """
Compute the wavenumbers of the drumhead.
Parameters:
----------
n_max_modes: int
The number of angular modes.
m_max_modes: int
The number of radial modes.
radius: float
The radius of the drumhead.
Returns:
-------
wavenumbers: np.ndarray
The wavenumbers for the drumhead.
"""
# Bessel function roots for different orders (n) and radial modes (m)
= np.zeros((n_max_modes, m_max_modes))
wavenumbers for n in range(n_max_modes):
= jn_zeros(n, m_max_modes)
bessel_roots = bessel_roots / radius
wavenumbers[n, :] return wavenumbers
def drumhead_eigenvalues(
# (n_max_modes, m_max_modes)
wavenumbers: np.ndarray,
):"""
Compute the eigenvalues of the drumhead.
The eigenvalues of the drumhead are given by the square of the wavenumbers.
Parameters:
----------
wavenumbers: np.ndarray
The wavenumbers for the drumhead.
squared: bool
If True, return the squared eigenvalues.
Returns:
-------
eigenvalues: np.ndarray
The eigenvalues of the drumhead.
"""
return wavenumbers**2
def drumhead_eigenfunctions(
# (n_max_modes, m_max_modes)
wavenumbers: np.ndarray, # (n_gridpoints_r)
r: np.ndarray, # (n_gridpoints_theta)
theta: np.ndarray, -> np.ndarray:
) """
Compute the modes of the drumhead.
The modes of the drumhead are given by the Bessel function times the sine/cosine of the angle.
Parameters:
----------
wavenumbers: np.ndarray
The wavenumbers for the drumhead.
r: np.ndarray
Radial grid points.
theta: np.ndarray
Angular grid points.
Returns:
-------
modes: np.ndarray
The eigenfunctions for the drumhead.
"""
= wavenumbers.shape
n_max_modes, m_max_modes = len(r)
n_gridpoints_r = len(theta)
n_gridpoints_theta
# Initialize the modes array
= np.zeros((n_max_modes, m_max_modes, n_gridpoints_r, n_gridpoints_theta))
modes = np.zeros_like(modes)
inverse_modes
= np.zeros((n_max_modes, m_max_modes))
squared_norms
= np.meshgrid(r, theta, indexing="ij")
r_grid, theta_grid
# Compute the modes
for n in range(n_max_modes):
for m in range(m_max_modes):
= np.cos(n * theta_grid) * jv(n, wavenumbers[n, m] * r_grid)
modes[n, m]
# calculate the squared norm
= modes[n, m] ** 2 * r_grid
integrand = simpson(integrand, x=theta, axis=1)
integral_theta = simpson(integral_theta, x=r)
squared_norms[n, m]
# get the normalised modes
= modes[n, m] / squared_norms[n, m]
inverse_modes[n, m]
return modes, inverse_modes, squared_norms
:::
::: {#cell-10 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def dblintegral(integrand, x, y, method="simpson"):
"""
Compute the double integral of a function K over the domain x and y.
"""
if method == "simpson":
= simpson(integrand, x=y, axis=1)
integral_y return simpson(integral_y, x=x)
elif method == "trapezoid":
= trapezoid(integrand, x=y, axis=1)
integral_y return trapezoid(integral_y, x=x)
else:
raise ValueError("Method not supported")
:::
# Example usage
= 25
n_max_modes = 25
m_max_modes = 1.0
radius = 100
n_gridpoints_r = 100
n_gridpoints_theta
= drumhead_wavenumbers(n_max_modes, m_max_modes, radius)
wavenumbers = drumhead_eigenvalues(wavenumbers)
eigenvalues = np.linspace(0, radius, n_gridpoints_r)
r = np.linspace(0, 2 * np.pi, n_gridpoints_theta)
theta = drumhead_eigenfunctions(wavenumbers, r, theta)
K_fwd, K_inv, K_N
assert K_inv.shape == (
n_max_modes,
m_max_modes,
n_gridpoints_r,
n_gridpoints_theta,
)
assert K_fwd.shape == (
n_max_modes,
m_max_modes,
n_gridpoints_r,
n_gridpoints_theta, )
\[ K_{n,m}(r, \varphi) = \cos (n \varphi) J_n\left(\mu_{n, m} \frac{r}{R}\right) \]
where \(J_n\) is the Bessel function of the first kind of order \(n\), and \(\mu_{n, m}\) is the \(m\)-th root of the \(n\)-th order Bessel function of the first kind.
= 6
n_max_modes_x = 6
n_max_modes_y = 20
n_gridpoints_x = 20
n_gridpoints_y = 1.08
length_x = 1.08
length_y = np.linspace(0, length_x, n_gridpoints_x)
grid_x = np.linspace(0, length_y, n_gridpoints_y)
grid_y
# slow version
# Define the range for n and m
= np.arange(1, n_max_modes_x + 1)
n_values = np.arange(1, n_max_modes_y + 1)
m_values
# Define the range for x and y
= np.linspace(0, length_x, n_gridpoints_x)
x_values = np.linspace(0, length_y, n_gridpoints_y)
y_values
# Initialize the 4D array to store the results
= np.zeros((len(n_values), len(m_values), len(x_values), len(y_values)))
K = np.zeros((len(n_values), len(m_values)))
Lambda # Compute the values
for i, n in enumerate(n_values):
for j, m in enumerate(m_values):
= np.pi**2 * ((n / length_x) ** 2 + (m / length_y) ** 2)
Lambda[i, j] for k, x in enumerate(x_values):
for l, y in enumerate(y_values):
= np.sin(n * np.pi * x / length_x) * np.sin(
K[i, j, k, l] * np.pi * y / length_y
m )
::: {#cell-14 .cell 0=‘t’ 1=‘e’ 2=‘s’ 3=‘t’}
= plate_wavenumbers(n_max_modes_x, n_max_modes_y, length_x, length_y)
wnx, wny assert np.allclose(plate_eigenfunctions(wnx, wny, grid_x, grid_y), K)
assert np.allclose(plate_eigenvalues(wnx, wny), Lambda)
:::
::: {#cell-15 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def forward_STL(
# (n_modes, n_gridpoints)
K: np.ndarray, # (n_gridpoints, n_samples) or (n_gridpoints,)
u: np.ndarray, float, # grid spacing
dx: -> np.ndarray:
) """
Compute the forward STL transform. The integration is done using the trapezoidal rule.
Parameters
----------
K: np.ndarray
The sampled eigenfunctions of the string. Shape (n_modes, n_gridpoints)
u: np.ndarray
The input signal. Shape (n_gridpoints, n_samples) or (n_gridpoints,)
dx: float
The grid spacing of the sampled string.
Returns
-------
np.ndarray
The transformed signal. Shape (n_modes, n_samples) or (n_modes,)
"""
if u.ndim == 1:
= u[:, None]
u = dx * einsum("m n, n s -> m s", K, u)
transformed_signal return transformed_signal if u.shape[1] > 1 else transformed_signal[:, 0]
def inverse_STL(
# (n_modes, n_gridpoints)
K: np.ndarray, # (n_modes, n_samples) or (n_modes,)
u_bar: np.ndarray, float, # length of the string
length: -> np.ndarray:
) """
Compute the inverse STL transform using the formula of Rabenstein et al. (2000).
Parameters:
-----------
K: np.ndarray
The sampled eigenfunctions of the string. Shape (n_modes, n_gridpoints)
u_bar: np.ndarray
The transformed signal. Shape (n_modes, n_samples) or (n_modes,)
length: float
The length of the string.
Returns:
--------
np.ndarray
The reconstructed signal. Shape (n_gridpoints, n_samples) or (n_gridpoints,)
"""
if u_bar.ndim == 1:
= u_bar[:, None] # Convert to (n_modes, 1) if input is (n_modes,)
u_bar
= length / 2.0
N = einsum("m n, m s -> n s", K, u_bar) / N
reconstructed_signal return reconstructed_signal if u_bar.shape[1] > 1 else reconstructed_signal[:, 0]
:::
::: {#cell-16 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def forward_STL_2d(
# (n_modes_x, n_modes_y, n_gridpoints_x, n_gridpoints_y)
K: np.ndarray, # (n_gridpoints_x, n_gridpoints_y, n_samples) or (n_gridpoints_x, n_gridpoints_y)
u: np.ndarray, float, # grid spacing
x: float, # grid spacing
y: bool = False,
use_simpson: -> np.ndarray:
) """
Compute the forward STL transform. The integration is done using the trapezoidal rule.
Parameters:
-----------
K: np.ndarray
The sampled eigenfunctions of the string. Shape (n_modes, n_gridpoints)
u: np.ndarray
The input signal. Shape (n_gridpoints, n_samples) or (n_gridpoints,)
dx: float
The grid spacing of the sampled string.
Returns:
--------
np.ndarray
The transformed signal. Shape (n_modes, n_samples) or (n_modes,)
"""
if u.ndim == 2:
= u[..., None]
u
if use_simpson:
= K.shape
n_modes_x, n_modes_y, _, _ = u.shape
_, _, n_samples
# print(u.shape)
= np.zeros((n_modes_x, n_modes_y, n_samples))
transformed_signal
for mode_x in range(n_modes_x):
for mode_y in range(n_modes_y):
for sample in range(n_samples):
# Perform 2D Simpson's integration
= K[mode_x, mode_y] * u[:, :, sample]
uu
# integral_x = simpson([simpson(uu_y, dx=dy) for uu_y in uu], dx=dx)
= dblintegral(
transformed_signal[mode_x, mode_y, sample]
uu,
x,
y,="simpson",
method
)# else use trapezoidal rule
else:
= x[1] - x[0]
dx = y[1] - y[0]
dy = dx * dy * einsum("m n x y, x y s -> m n s", K, u)
transformed_signal
return (
transformed_signal.squeeze()if transformed_signal.shape[-1] == 1
else transformed_signal
)
def inverse_STL_2d(
# (n_modes_x, n_modes_y, n_gridpoints_x, n_gridpoints_y)
K: np.ndarray, # (n_modes_x, n_modes_y, n_samples) or (n_modes_x, n_modes_y)
u_bar: np.ndarray, float, # length in x
l1: float, # length in y
l2: -> np.ndarray:
) """
Compute the inverse STL transform using the formula of Rabenstein et al. (2000).
Parameters:
-----------
K: np.ndarray
The sampled eigenfunctions of the string. Shape (n_modes, n_gridpoints)
u_bar: np.ndarray
The transformed signal. Shape (n_modes, n_samples) or (n_modes,)
length: float
The length of the string.
Returns:
--------
np.ndarray
The reconstructed signal. Shape (n_gridpoints, n_samples) or (n_gridpoints,)
"""
if u_bar.ndim == 2:
= u_bar[
u_bar None
..., # Convert to (n_modes_x, n_modes_y, 1) if input is (n_modes_x, n_modes_y)
]
= 0.25 * (l1 * l2)
N = einsum("m n x y, m n s -> x y s", K, u_bar) / N
reconstructed_signal return (
if u_bar.shape[-1] == 1 else reconstructed_signal
reconstructed_signal.squeeze() )
:::
= 1.08
length_x = 0.8
length_y = 25
n_max_modes_x = 25
n_max_modes_y = 100
n_gridpoints_x = 100
n_gridpoints_y
= np.linspace(0, length_x, n_gridpoints_x)
x = np.linspace(0, length_y, n_gridpoints_y)
y
= plate_wavenumbers(
wnx, wny
n_max_modes_x,
n_max_modes_y,
length_x,
length_y,
)= plate_eigenfunctions(wnx, wny, x, y)
K
= 0.5 * K[2, 2] + 0.5 * K[3, 3]
g
= forward_STL_2d(K, g, x, y, use_simpson=True)
bar_g = inverse_STL_2d(K, bar_g, length_x, length_y)
g_reconstructed
assert np.allclose(g, g_reconstructed, atol=1e-2)
= plt.subplots(1, 2, figsize=(10, 5))
fig, ax
0].imshow(g, origin="lower", aspect="auto")
ax[0].set_title("Original excitation")
ax[1].imshow(g_reconstructed, origin="lower", aspect="auto")
ax[1].set_title("Reconstructed excitation") ax[
::: {#cell-18 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def forward_STL_drumhead(
# (n_modes_r, n_modes_theta, n_gridpoints_r, n_gridpoints_theta)
K: np.ndarray, # (n_gridpoints_x, n_gridpoints_y, n_samples) or (n_gridpoints_x, n_gridpoints_y)
u: np.ndarray, # radial grid
r: np.ndarray, # angular grid
theta: np.ndarray, bool = False,
use_simpson: -> np.ndarray:
) """
Compute the forward STL transform. The integration is done using the trapezoidal rule or Simpson's rule.
Parameters:
-----------
K: np.ndarray
The sampled eigenfunctions of the string. Shape (n_modes_r, n_modes_theta, n_gridpoints_r, n_gridpoints_theta)
u: np.ndarray
The input signal. Shape (n_gridpoints_r, n_gridpoints_theta, n_samples) or (n_gridpoints_r, n_gridpoints_theta)
dr: float
The grid spacing of the sampled membrane.
dtheta: float
The grid spacing of the sampled membrane.
Returns:
--------
np.ndarray
The transformed signal. Shape (n_modes_r, n_modes_theta, n_samples) or (n_modes_r, n_modes_theta)
"""
if u.ndim == 2:
= u[..., None]
u
= np.meshgrid(r, theta, indexing="ij")
r_grid, _
if use_simpson:
= K.shape
max_n, max_m, _, _ = u.shape
_, _, n_samples
= np.zeros((max_n, max_m, n_samples))
transformed_signal
for n in range(max_n):
for m in range(max_m):
for sample in range(n_samples):
= (
integrand * K[n, m] * r_grid
u[..., sample] # notice the r_grid (Jacobian determinant)
) = dblintegral(
transformed_signal[n, m] =r, y=theta, method="trapezoid"
integrand, x
)
else:
# integrand has shape (n_modes_r, n_modes_theta, n_gridpoints_r, n_gridpoints_theta)
= u[..., None].transpose(2, 3, 0, 1) * K * r_grid
integrand
= trapezoid(integrand, x=theta, axis=-1)
integral_y = trapezoid(integral_y, x=r, axis=-1)
transformed_signal
return (
transformed_signal.squeeze()if transformed_signal.shape[-1] == 1
else transformed_signal
)
def inverse_STL_drumhead(
# (n_modes_x, n_modes_y, n_gridpoints_x, n_gridpoints_y)
K_inv: np.ndarray, # (n_modes_x, n_modes_y, n_samples) or (n_modes_x, n_modes_y)
u_bar: np.ndarray, -> np.ndarray:
) """
Compute the inverse STL transform using the formula of Rabenstein et al. (2000).
Parameters:
-----------
K: np.ndarray
The sampled eigenfunctions of the string. Shape (n_modes, n_gridpoints)
u_bar: np.ndarray
The transformed signal. Shape (n_modes, n_samples) or (n_modes,)
length: float
The length of the string.
Returns:
--------
np.ndarray
The reconstructed signal. Shape (n_gridpoints, n_samples) or (n_gridpoints,)
"""
if u_bar.ndim == 2:
= u_bar[
u_bar None
..., # Convert to (n_modes_x, n_modes_y, 1) if input is (n_modes_x, n_modes_y)
]
= einsum("n m x y, n m s -> x y s", K_inv, u_bar)
reconstructed_signal return (
if u_bar.shape[-1] == 1 else reconstructed_signal
reconstructed_signal.squeeze() )
:::
# Example usage
= 10
n_max_modes = 10
m_max_modes = 1.0
radius = 100
n_gridpoints_r = 100
n_gridpoints_theta
= drumhead_wavenumbers(n_max_modes, m_max_modes, radius)
wavenumbers = drumhead_eigenvalues(wavenumbers)
eigenvalues = np.linspace(0, radius, n_gridpoints_r)
r = np.linspace(0, 2 * np.pi, n_gridpoints_theta)
theta = drumhead_eigenfunctions(wavenumbers, r, theta)
K_fwd, K_inv, K_N
assert np.allclose(
K_fwd.shape, (n_max_modes, m_max_modes, n_gridpoints_r, n_gridpoints_theta)# Should be (10, 10, 100, 100)
) assert np.allclose(
K_inv.shape, (n_max_modes, m_max_modes, n_gridpoints_r, n_gridpoints_theta)# Should be (10, 10, 100, 100)
)
# Create an example g array to test the transforms
= K_fwd[3, 3]
g
= forward_STL_drumhead(K_fwd, g, r, theta, use_simpson=False)
bar_g = inverse_STL_drumhead(K_inv, bar_g)
g_reconstructed
# Verify if g can be reconstructed
assert np.allclose(g, g_reconstructed, atol=1e-2)
print(g.min(), g.max())
print(g_reconstructed.min(), g_reconstructed.max())
# Plot using pcolormesh
= plt.subplots(
fig, ax 1,
2,
={"projection": "polar"},
subplot_kw=(10, 5),
figsize
)= ax[0].pcolormesh(theta, r, g, shading="auto", cmap="viridis")
c = ax[1].pcolormesh(theta, r, g_reconstructed, shading="auto", cmap="viridis") c
::: {#cell-20 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def eigenvalues_from_pde(
strpars: StringParameters,
wavenumbers: np.ndarray,-> np.ndarray:
) """
Compute the positive imaginary side of the eigenvalues of the continuous-time system from the PDE parameters.
Parameters
----------
d1 : float
The linear damping coefficient.
d3 : float
The cubic damping coefficient.
rho : float
The density of the string.
A : float
The cross-sectional area of the string.
E : float
The Young's modulus of the string.
I : float
The second moment of area of the string.
wavenumbers : np.ndarray
The wavenumbers of the modes.
Returns
-------
np.ndarray
The eigenvalues of the continuous-time system.
"""
= strpars.d1 + strpars.d3 * wavenumbers**2
sigma_mu = sigma_mu / (2 * strpars.area_density)
sigma_mu
= strpars.E * strpars.I * wavenumbers**4 + strpars.Ts0 * wavenumbers**2
beta_mu = (beta_mu / strpars.area_density) - sigma_mu**2
omega_mu_squared
= np.sqrt(omega_mu_squared) # imaginary-part frequency
omega_mu
= -sigma_mu + 1j * omega_mu # continuous time eigenvalues
gt_eigenvalues return gt_eigenvalues
:::
from IPython.display import Audio
= 50
n_max_modes = 44100
sr = 1 / sr
dt = 1.0
final_time = int(final_time / dt)
n_samples = StringParameters()
p_params = string_eigenvalues_sqrt(n_max_modes, p_params.length)
wn = eigenvalues_from_pde(p_params, wn)
eigvals
print(eigvals.imag.min() / (2 * np.pi))
= np.exp(eigvals * dt)
eigvals_d
= np.vander(eigvals_d, n_samples, increasing=True).real
states sum(0), rate=sr)) display(Audio(states.
::: {#cell-23 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def eigenvalues_from_plate_pde(
platepars: PlateParameters,# wavenumbers x (n_max_modes_x,)
wnx: np.ndarray, # wavenumbers y (n_max_modes_y,)
wny: np.ndarray, -> np.ndarray:
) """
Compute the positive imaginary side of the eigenvalues of the continuous-time system from the PDE parameters of the rectangular plate. From 5.96 of Digital Sound Synthesis using the FTM, and Eq. 8 of TENSION MODULATED NONLINEAR 2D MODELS FOR DIGITAL SOUND SYNTHESIS WITH THE FUNCTIONAL TRANSFORMATION METHOD.
Parameters
----------
d1 : float
The frequency independent damping coefficient.
d3 : float
The frequency dependent damping coefficient.
rho : float
The density of the string.
A : float
The cross-sectional area of the string.
E : float
The Young's modulus of the string.
I : float
The second moment of area of the string.
wavenumbers : np.ndarray
The wavenumbers of the modes.
Returns
-------
np.ndarray
The eigenvalues of the continuous-time system.
"""
# The term that goes with the fourth spatial derivative
= platepars.E * platepars.h**3 / (12 * (1 - platepars.nu**2))
S # S = platepars.E * platepars.I / (1 - platepars.nu)
# Lambda_squared = plate_eigenvalues(wnx, wny, squared=True)
= plate_eigenvalues(wnx, wny)
Lambda
= platepars.d0 + platepars.d2 * Lambda # real-part decay
sigma_mu = sigma_mu / (2 * platepars.surface_density)
sigma_mu
= (
beta_mu * Lambda + S * Lambda**2
platepars.Tm # this gives a correct result but it is not the same as Lambda_squared
) = beta_mu / platepars.surface_density - sigma_mu**2
omega_mu_squared
= np.sqrt(omega_mu_squared) # imaginary-part frequency
omega_mu
= -sigma_mu + 1j * omega_mu # continuous time eigenvalues
gt_eigenvalues return gt_eigenvalues
:::
= 20
n_max_modes_x = 20
n_max_modes_y = 44100
sr = 1 / sr
dt = 6.0
final_time = int(final_time / dt)
n_samples
= PlateParameters()
p_params = plate_wavenumbers(n_max_modes_x, n_max_modes_y, p_params.l1, p_params.l2)
wnx, wny = eigenvalues_from_plate_pde(p_params, wnx, wny).reshape(-1)
eigvals
print(eigvals.imag.min() / (2 * np.pi))
= np.exp(eigvals * dt)
eigvals_d
= np.vander(eigvals_d, n_samples, increasing=True).real
states
sum(0), rate=sr)) display(Audio(states.
::: {#cell-25 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def eigenvalues_from_drumhead_pde(
drumhead_pars: CircularDrumHeadParameters,# (n_max_modes_r, m_max_modes_theta)
Lambda_nm: np.ndarray, -> np.ndarray:
) """
Compute the positive imaginary side of the eigenvalues of the continuous-time system from the PDE parameters of the rectangular plate. From 5.96 of Digital Sound Synthesis using the FTM, and Eq. 8 of TENSION MODULATED NONLINEAR 2D MODELS FOR DIGITAL SOUND SYNTHESIS WITH THE FUNCTIONAL TRANSFORMATION METHOD.
Parameters
----------
d1 : float
The frequency independent damping coefficient.
d3 : float
The frequency dependent damping coefficient.
rho : float
The density of the string.
A : float
The cross-sectional area of the string.
E : float
The Young's modulus of the string.
I : float
The second moment of area of the string.
wavenumbers : np.ndarray
The wavenumbers of the modes.
Returns
-------
np.ndarray
The eigenvalues of the continuous-time system.
"""
# The term that goes with the fourth spatial derivative
= drumhead_pars.E * drumhead_pars.h**3 / (12 * (1 - drumhead_pars.nu**2))
D
= drumhead_pars.d0 + drumhead_pars.d2 * Lambda_nm # real-part decay
sigma_mu = sigma_mu / (2 * drumhead_pars.surface_density)
sigma_mu
= D * Lambda_nm**2 + drumhead_pars.Tm * Lambda_nm
omega_mu_squared = omega_mu_squared / drumhead_pars.surface_density - sigma_mu**2
omega_mu_squared
= np.sqrt(omega_mu_squared) # imaginary-part frequency
omega_mu
= -sigma_mu + 1j * omega_mu # continuous time eigenvalues
gt_eigenvalues return gt_eigenvalues
:::
= 20
n_max_modes_x = 20
n_max_modes_y = 44100
sr = 1 / sr
dt = 2.0
final_time = int(final_time / dt)
n_samples = CircularDrumHeadParameters.avanzini()
p_params
= drumhead_eigenvalues(
v
drumhead_wavenumbers(n_max_modes_x, n_max_modes_y, p_params.r0)
)= eigenvalues_from_drumhead_pde(p_params, v).reshape(-1)
eigvals
# import matplotlib.pyplot as plt
print(eigvals.imag.min() / (2 * np.pi))
print(eigvals.imag.max() / (2 * np.pi))
# plt.plot(eigvals)
= np.exp(eigvals * dt)
eigvals_d = np.vander(eigvals_d, n_samples, increasing=True).real
states
# fig, ax = plt.subplots(1, 1, figsize=(10, 5))
# ax.plot(states[0])
0], rate=sr)) display(Audio(states[
::: {#cell-27 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def sample_parallel_tf(
# (n_modes,)
num: np.ndarray, # (n_modes,)
den: np.ndarray, float,
dt: str = "impulse",
method:
):"""
Sample a parallel transfer function using the impulse invariant method.
Parameters
----------
num : np.ndarray
The numerator of the transfer function.
den : np.ndarray
The denominator of the transfer function.
Returns
-------
num_d : np.ndarray
The numerator of the discrete-time transfer function.
den_d : np.ndarray
The denominator of the discrete-time transfer function.
"""
def sample(n, d):
= cont2discrete((n, d), dt, method=method)
b, a, _ return b.flatten(), a.flatten()
= zip(*[sample(n, d) for n, d in zip(num, den)])
b, a
return np.array(b), np.array(a)
:::
::: {#cell-28 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def tf_excitation_continuous(
eigenvalues: np.ndarray,float, # surface or area density
density: -> tuple[np.ndarray, np.ndarray]:
) """
Compute the continuous excitation transfer function.
Parameters
----------
eigenvalues : np.ndarray
The eigenvalues of the system.
density: float
The surface density of a membrane (rho * h) or area density of a string (rho * A)
dt : float
The time step size.
Returns
-------
np.ndarray
The numerator of the discrete-time transfer function.
np.ndarray
The denominator of the discrete-time transfer function
"""
= -eigenvalues.real
sigma_mu = eigenvalues.imag
omega_mu
= sigma_mu * 2
a1 = sigma_mu**2 + omega_mu**2
a2
= np.ones_like(sigma_mu)
ones = ones / density
b = np.stack([ones, a1, a2], axis=-1)
a return b, a
def tf_excitation_discrete(
eigenvalues: np.ndarray,float, # surface or area density
density: float, # time step
dt: -> tuple[np.ndarray, np.ndarray]:
) """
Compute the discrete-time excitation transfer function of a system.
Parameters
----------
eigenvalues : np.ndarray
The eigenvalues of the system.
density: float
The surface density of a membrane (rho * h) or area density of a string (rho * A)
dt : float
The time step size.
Returns
-------
np.ndarray
The numerator of the discrete-time transfer function.
np.ndarray
The denominator of the discrete-time transfer function
"""
= tf_excitation_continuous(eigenvalues, density)
b, a
# Discretize the system
= sample_parallel_tf(b, a, dt)
tf_d
return tf_d
:::
::: {#cell-29 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
def tf_initial_conditions_continuous(
eigenvalues: np.ndarray,-> tuple[np.ndarray, np.ndarray]:
) """
Compute the continuos "initial-conditions" transfer function from the eigenvalues of the system.
Parameters
----------
eigenvalues : np.ndarray
The eigenvalues of the system.
density: float
The surface density of a membrane (rho * h) or area density of a string (rho * A)
dt : float
The time step size.
Returns
-------
np.ndarray
The numerator of the discrete-time transfer function.
np.ndarray
The denominator of the discrete-time transfer function
"""
= -eigenvalues.real
sigma_mu = eigenvalues.imag
omega_mu
= sigma_mu * 2
a1 = sigma_mu**2 + omega_mu**2
a2
= np.ones_like(sigma_mu)
ones = a1
b1
= np.stack([ones, b1], axis=-1)
b = np.stack([ones, a1, a2], axis=-1)
a
return (b, a)
def tf_initial_conditions_discrete(
eigenvalues: np.ndarray,float, # time step
dt: -> tuple[np.ndarray, np.ndarray]:
) """
Compute the discrete-time initial conditions transfer function of a system.
Parameters
----------
eigenvalues : np.ndarray
The eigenvalues of the system.
density: float
The surface density of a membrane (rho * h) or area density of a string (rho * A)
dt : float
The time step size.
Returns
-------
np.ndarray
The numerator of the discrete-time transfer function.
np.ndarray
The denominator of the discrete-time transfer function
"""
= tf_initial_conditions_continuous(eigenvalues)
b, a
# Discretize the system
= sample_parallel_tf(b, a, dt)
tf_d
return tf_d
:::
::: {#cell-30 .cell 0=‘t’ 1=‘e’ 2=‘s’ 3=‘t’}
= tf_excitation_discrete(eigvals, p_params.surface_density, dt)
b, a = tf_initial_conditions_discrete(eigvals, dt)
b_ic, a_ic
# manual discretization
= np.exp(eigvals * dt)
eigenvalues_d
# for the excitation tf
= (
b1 * dt)
np.exp(eigvals.real * np.sin(eigvals.imag * dt)
/ eigvals.imag
/ p_params.surface_density
)
# for the initial conditions tf
# here we ignore initial velocity
= np.exp(eigvals.real * dt)
r = r * np.sin(eigvals.imag * dt) / eigvals.imag * -eigvals.real - r * np.cos(
b1_ic * dt
eigvals.imag
)
= -2 * np.exp(eigvals.real * dt) * np.cos(eigvals.imag * dt)
a1 = np.exp(2 * eigvals.real * dt)
a2 = np.stack([np.zeros_like(b1), b1, np.zeros_like(b1)], axis=-1) * dt
b_manual = np.stack([np.ones_like(a1), a1, a2], axis=-1)
a_manual
= np.stack([np.ones_like(b1_ic), b1_ic], axis=-1) * dt
b_ic_manual
print(b[0])
print(b_manual[0])
print(b_ic[0] * sr)
print(b_ic_manual[0] * sr)
assert np.allclose(b[:, 1], b_manual[:, 1])
assert np.allclose(a, a_manual)
assert np.allclose(b_ic[:, :2], b_ic_manual[:, :2])
:::