import matplotlib.pyplot as plt
Funtional Transformation Method Utilities
This notebook contains a set of utility functions for the functional transformation method.
CircularDrumHeadParameters
CircularDrumHeadParameters (h:float=0.00019, r0:float=0.328, I:float=5.7e-13, rho:float=1380.0, E:int=3500000000.0, nu:float=0.35, d0:float=0.14, d2:float=0.32, Tm:float=3990, f0:float=143.95)
Kettle drum head, from Digital sound synthesis of string instruments with the functional transformation method Table 5.2.
PlateParameters
PlateParameters (h:float=0.0012, l1:float=1.08, l2:float=1.08, I:float=1.44e-10, rho:float=7800.0, E:int=220000000000.0, nu:float=0.29, d0:float=0.032, d2:float=0.0013, Tm:float=2010, f_1:float=10.27)
From Digital sound synthesis of string instruments with the functional transformation method Table 5.2.
StringParameters
StringParameters (A:float=5.188e-07, I:float=1.41e-13, rho:float=1140, E:int=5400000000.0, d1:float=8e-05, d3:float=1.4e-05, Ts0:float=60.97, length:float=0.65)
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.
string_eigenfunctions
string_eigenfunctions (wavenumbers:numpy.ndarray, grid:numpy.ndarray)
*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.*
string_eigenvalues_sqrt
string_eigenvalues_sqrt (n_max_modes:int, length:float)
*Compute the square root of the eigenvalues of the string (in this case the wavenumbers). The eigenvalues of the string are given by:
\[ \lambda_{\mu} = (rac{\mu \pi}{L}) \]
where \(\mu\) is the mode number and \(L\) is the length of the string.*
plate_eigenfunctions
plate_eigenfunctions (wavenumbers_x:numpy.ndarray, wavenumbers_y:numpy.ndarray, x:numpy.ndarray, y:numpy.ndarray)
*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.*
Type | Details | |
---|---|---|
wavenumbers_x | ndarray | (n_max_modes_x,) |
wavenumbers_y | ndarray | (n_max_modes_y,) |
x | ndarray | (n_gridpoints_x,) |
y | ndarray | (n_gridpoints_y,) |
Returns | ndarray |
plate_eigenvalues
plate_eigenvalues (wavenumbers_x:numpy.ndarray, wavenumbers_y:numpy.ndarray, squared:bool=False)
*Compute the eigenvalues of the plate. The eigenvalues of the plate are given by:
$$
\lambda_{\mu,
u} = $$
where $\mu$ and $
u$ are the mode numbers and \(L_1\) and \(L_2\) are the width and height of the plate.*
Type | Default | Details | |
---|---|---|---|
wavenumbers_x | ndarray | (n_max_modes_x,) | |
wavenumbers_y | ndarray | (n_max_modes_y,) | |
squared | bool | False |
/Users/diaz/mambaforge/envs/physmodjax/lib/python3.10/site-packages/fastcore/docscrape.py:230: UserWarning: potentially wrong underline length...
Parameters:
---------- in
Compute the wavenumbers of the plate.
...
else: warn(msg)
/Users/diaz/mambaforge/envs/physmodjax/lib/python3.10/site-packages/fastcore/docscrape.py:230: UserWarning: potentially wrong underline length...
Returns:
------- in
Compute the wavenumbers of the plate.
...
else: warn(msg)
plate_wavenumbers
plate_wavenumbers (n_max_modes_x:int, n_max_modes_y:int, l1:float, l2:float)
*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*
/Users/diaz/mambaforge/envs/physmodjax/lib/python3.10/site-packages/fastcore/docscrape.py:230: UserWarning: potentially wrong underline length...
Parameters:
---------- in
Compute the modes of the drumhead.
The modes of the drumhead are given by the Bessel function times the sine/cosine of the angle....
else: warn(msg)
/Users/diaz/mambaforge/envs/physmodjax/lib/python3.10/site-packages/fastcore/docscrape.py:230: UserWarning: potentially wrong underline length...
Returns:
------- in
Compute the modes of the drumhead.
The modes of the drumhead are given by the Bessel function times the sine/cosine of the angle....
else: warn(msg)
drumhead_eigenfunctions
drumhead_eigenfunctions (wavenumbers:numpy.ndarray, r:numpy.ndarray, theta:numpy.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.*
Type | Details | |
---|---|---|
wavenumbers | ndarray | (n_max_modes, m_max_modes) |
r | ndarray | (n_gridpoints_r) |
theta | ndarray | (n_gridpoints_theta) |
Returns | ndarray |
/Users/diaz/mambaforge/envs/physmodjax/lib/python3.10/site-packages/fastcore/docscrape.py:230: UserWarning: potentially wrong underline length...
Parameters:
---------- in
Compute the eigenvalues of the drumhead.
The eigenvalues of the drumhead are given by the square of the wavenumbers....
else: warn(msg)
/Users/diaz/mambaforge/envs/physmodjax/lib/python3.10/site-packages/fastcore/docscrape.py:230: UserWarning: potentially wrong underline length...
Returns:
------- in
Compute the eigenvalues of the drumhead.
The eigenvalues of the drumhead are given by the square of the wavenumbers....
else: warn(msg)
drumhead_eigenvalues
drumhead_eigenvalues (wavenumbers:numpy.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.*
Type | Details | |
---|---|---|
wavenumbers | ndarray | (n_max_modes, m_max_modes) |
/Users/diaz/mambaforge/envs/physmodjax/lib/python3.10/site-packages/fastcore/docscrape.py:230: UserWarning: potentially wrong underline length...
Parameters:
---------- in
Compute the wavenumbers of the drumhead.
...
else: warn(msg)
/Users/diaz/mambaforge/envs/physmodjax/lib/python3.10/site-packages/fastcore/docscrape.py:230: UserWarning: potentially wrong underline length...
Returns:
------- in
Compute the wavenumbers of the drumhead.
...
else: warn(msg)
drumhead_wavenumbers
drumhead_wavenumbers (n_max_modes:int, m_max_modes:int, radius:float)
*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.*
dblintegral
dblintegral (integrand, x, y, method='simpson')
Compute the double integral of a function K over the domain x and y.
# 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(m * np.pi * y / length_y) K[i, j, k, l]
= 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)
/Users/diaz/mambaforge/envs/physmodjax/lib/python3.10/site-packages/fastcore/docscrape.py:230: UserWarning: Unknown section Parameters:
else: warn(msg)
/Users/diaz/mambaforge/envs/physmodjax/lib/python3.10/site-packages/fastcore/docscrape.py:230: UserWarning: Unknown section Returns:
else: warn(msg)
inverse_STL
inverse_STL (K:numpy.ndarray, u_bar:numpy.ndarray, length:float)
Compute the inverse STL transform using the formula of Rabenstein et al. (2000).
Type | Details | |
---|---|---|
K | ndarray | (n_modes, n_gridpoints) |
u_bar | ndarray | (n_modes, n_samples) or (n_modes,) |
length | float | length of the string |
Returns | ndarray |
forward_STL
forward_STL (K:numpy.ndarray, u:numpy.ndarray, dx:float)
Compute the forward STL transform. The integration is done using the trapezoidal rule.
Type | Details | |
---|---|---|
K | ndarray | (n_modes, n_gridpoints) |
u | ndarray | (n_gridpoints, n_samples) or (n_gridpoints,) |
dx | float | grid spacing |
Returns | ndarray |
inverse_STL_2d
inverse_STL_2d (K:numpy.ndarray, u_bar:numpy.ndarray, l1:float, l2:float)
Compute the inverse STL transform using the formula of Rabenstein et al. (2000).
Type | Details | |
---|---|---|
K | ndarray | (n_modes_x, n_modes_y, n_gridpoints_x, n_gridpoints_y) |
u_bar | ndarray | (n_modes_x, n_modes_y, n_samples) or (n_modes_x, n_modes_y) |
l1 | float | length in x |
l2 | float | length in y |
Returns | ndarray |
forward_STL_2d
forward_STL_2d (K:numpy.ndarray, u:numpy.ndarray, x:float, y:float, use_simpson:bool=False)
Compute the forward STL transform. The integration is done using the trapezoidal rule.
Type | Default | Details | |
---|---|---|---|
K | ndarray | (n_modes_x, n_modes_y, n_gridpoints_x, n_gridpoints_y) | |
u | ndarray | (n_gridpoints_x, n_gridpoints_y, n_samples) or (n_gridpoints_x, n_gridpoints_y) | |
x | float | grid spacing | |
y | float | grid spacing | |
use_simpson | bool | False | |
Returns | ndarray |
= 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[
Text(0.5, 1.0, 'Reconstructed excitation')
inverse_STL_drumhead
inverse_STL_drumhead (K_inv:numpy.ndarray, u_bar:numpy.ndarray)
Compute the inverse STL transform using the formula of Rabenstein et al. (2000).
Type | Details | |
---|---|---|
K_inv | ndarray | (n_modes_x, n_modes_y, n_gridpoints_x, n_gridpoints_y) |
u_bar | ndarray | (n_modes_x, n_modes_y, n_samples) or (n_modes_x, n_modes_y) |
Returns | ndarray |
forward_STL_drumhead
forward_STL_drumhead (K:numpy.ndarray, u:numpy.ndarray, r:numpy.ndarray, theta:numpy.ndarray, use_simpson:bool=False)
Compute the forward STL transform. The integration is done using the trapezoidal rule or Simpson’s rule.
Type | Default | Details | |
---|---|---|---|
K | ndarray | (n_modes_r, n_modes_theta, n_gridpoints_r, n_gridpoints_theta) | |
u | ndarray | (n_gridpoints_x, n_gridpoints_y, n_samples) or (n_gridpoints_x, n_gridpoints_y) | |
r | ndarray | radial grid | |
theta | ndarray | angular grid | |
use_simpson | bool | False | |
Returns | ndarray |
# 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
-0.4320502477222054 0.43401550352665574
-0.4320519880367027 0.4340172517572766
eigenvalues_from_pde
eigenvalues_from_pde (strpars:__main__.StringParameters, wavenumbers:numpy.ndarray)
Compute the positive imaginary side of the eigenvalues of the continuous-time system from the PDE parameters.
Type | Details | |
---|---|---|
strpars | StringParameters | |
wavenumbers | ndarray | The wavenumbers of the modes. |
Returns | ndarray | The eigenvalues of the continuous-time system. |
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 = np.arange(1, n_max_modes + 1)
mu = np.pi * mu / 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.
247.0158068963442
eigenvalues_from_plate_pde
eigenvalues_from_plate_pde (platepars:__main__.PlateParameters, wnx:numpy.ndarray, wny:numpy.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.
Type | Details | |
---|---|---|
platepars | PlateParameters | |
wnx | ndarray | wavenumbers x (n_max_modes_x,) |
wny | ndarray | wavenumbers y (n_max_modes_y,) |
Returns | ndarray | The eigenvalues of the continuous-time system. |
= 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.
10.456826975580865
eigenvalues_from_drumhead_pde
eigenvalues_from_drumhead_pde (drumhead_pars:__main__.CircularDrumHeadPa rameters, Lambda_nm:numpy.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.
Type | Details | |
---|---|---|
drumhead_pars | CircularDrumHeadParameters | |
Lambda_nm | ndarray | (n_max_modes_r, m_max_modes_theta) |
Returns | ndarray | The eigenvalues of the continuous-time system. |
= 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(drumhead_wavenumbers(n_max_modes_x, n_max_modes_y, p_params.r0))
v = 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[
142.6549044721433
6141.459247890741
sample_parallel_tf
sample_parallel_tf (num:numpy.ndarray, den:numpy.ndarray, dt:float, method:str='impulse')
Sample a parallel transfer function using the impulse invariant method.
Type | Default | Details | |
---|---|---|---|
num | ndarray | (n_modes,) | |
den | ndarray | (n_modes,) | |
dt | float | ||
method | str | impulse | |
Returns | np.ndarray | The numerator of the discrete-time transfer function. |
tf_excitation_discrete
tf_excitation_discrete (eigenvalues:numpy.ndarray, density:float, dt:float)
Compute the discrete-time excitation transfer function of a system.
Type | Details | |
---|---|---|
eigenvalues | ndarray | The eigenvalues of the system. |
density | float | surface or area density |
dt | float | time step |
Returns | Tuple | The numerator of the discrete-time transfer function. |
tf_excitation_continuous
tf_excitation_continuous (eigenvalues:numpy.ndarray, density:float)
Compute the continuous excitation transfer function.
Type | Details | |
---|---|---|
eigenvalues | ndarray | The eigenvalues of the system. |
density | float | surface or area density |
Returns | Tuple | The numerator of the discrete-time transfer function. |
tf_initial_conditions_discrete
tf_initial_conditions_discrete (eigenvalues:numpy.ndarray, dt:float)
Compute the discrete-time initial conditions transfer function of a system.
Type | Details | |
---|---|---|
eigenvalues | ndarray | The eigenvalues of the system. |
dt | float | time step |
Returns | Tuple | The numerator of the discrete-time transfer function. |
tf_initial_conditions_continuous
tf_initial_conditions_continuous (eigenvalues:numpy.ndarray)
Compute the continuos “initial-conditions” transfer function from the eigenvalues of the system.
Type | Details | |
---|---|---|
eigenvalues | ndarray | The eigenvalues of the system. |
Returns | Tuple | The numerator of the discrete-time transfer function. |
= 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
= (
b1_ic * dt)
np.exp(eigvals.real * (
* dt) / eigvals.imag * -eigvals.real
np.sin(eigvals.imag - np.cos(eigvals.imag * dt)
)
)
= -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])
assert np.allclose(b[:, 1], b_manual[:, 1])
assert np.allclose(a, a_manual)
assert np.allclose(b_ic[:, :2], b_ic_manual[:, :2])
[0.00000000e+00 1.90416682e-09 1.11022302e-16]
[0.00000000e+00 1.90416703e-09 0.00000000e+00]