commit 785a780798483f8f086a47e7ce48d905b2612b7b
parent be7f9f1d8db68ad55ec89c80a49f3192d767d962
Author: minerva-jupiter <ryouturn@gmail.com>
Date: Wed, 17 Jun 2026 21:33:12 +0900
feat: implement physical bore resonance and reed coupling
- Introduce `Bore` struct for digital waveguide physical modeling, featuring non-linear dissipative boundary conditions
- Implement `Reed` struct for mass-spring-damper oscillator coupling
- Integrate `rustfft` for frequency domain analysis in CLI
- Update `Aerothesis` plugin to use full physical modeling synthesis instead of simple oscillation
- Add `OscillationType` and `BoundaryType` enums for user-configurable synthesis modes
- Update project metadata to reflect transition to an instrument plugin
Diffstat:
| M | Cargo.lock | | | 73 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| M | Cargo.toml | | | 1 | + |
| M | README.md | | | 221 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------- |
| M | src/lib.rs | | | 289 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--- |
| M | src/main.rs | | | 76 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------ |
5 files changed, 575 insertions(+), 85 deletions(-)
diff --git a/Cargo.lock b/Cargo.lock
@@ -22,6 +22,7 @@ name = "aerothesis"
version = "0.1.0"
dependencies = [
"nih_plug",
+ "rustfft",
"textplots",
]
@@ -70,6 +71,12 @@ dependencies = [
]
[[package]]
+name = "autocfg"
+version = "1.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53"
+
+[[package]]
name = "backtrace"
version = "0.3.76"
source = "registry+https://github.com/rust-lang/crates.io-index"
@@ -431,12 +438,39 @@ dependencies = [
]
[[package]]
+name = "num-complex"
+version = "0.4.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
name = "num-conv"
version = "0.2.2"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "521739c6d2bac4aa25192232afe6841231376b2b26d4d9fae5ecf8ca5772e441"
[[package]]
+name = "num-integer"
+version = "0.1.46"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-traits"
+version = "0.2.19"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841"
+dependencies = [
+ "autocfg",
+]
+
+[[package]]
name = "num_threads"
version = "0.1.7"
source = "registry+https://github.com/rust-lang/crates.io-index"
@@ -505,6 +539,15 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "439ee305def115ba05938db6eb1644ff94165c5ab5e9420d1c1bcedbba909391"
[[package]]
+name = "primal-check"
+version = "0.3.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dc0d895b311e3af9902528fbb8f928688abbd95872819320517cc24ca6b2bd08"
+dependencies = [
+ "num-integer",
+]
+
+[[package]]
name = "proc-macro2"
version = "1.0.106"
source = "registry+https://github.com/rust-lang/crates.io-index"
@@ -562,6 +605,20 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "b50b8869d9fc858ce7266cce0194bd74df58b9d0e3f6df3a9fc8eb470d95c09d"
[[package]]
+name = "rustfft"
+version = "6.4.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "21db5f9893e91f41798c88680037dba611ca6674703c1a18601b01a72c8adb89"
+dependencies = [
+ "num-complex",
+ "num-integer",
+ "num-traits",
+ "primal-check",
+ "strength_reduce",
+ "transpose",
+]
+
+[[package]]
name = "scopeguard"
version = "1.2.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
@@ -656,6 +713,12 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "8ed6a63f02c8539c91a8685a86f4099661ba3da017932f6ebbea6de3f0fa7c90"
[[package]]
+name = "strength_reduce"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "fe895eb47f22e2ddd4dabc02bce419d2e643c8e3b585c78158b349195bc24d82"
+
+[[package]]
name = "syn"
version = "1.0.109"
source = "registry+https://github.com/rust-lang/crates.io-index"
@@ -784,6 +847,16 @@ dependencies = [
]
[[package]]
+name = "transpose"
+version = "0.2.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1ad61aed86bc3faea4300c7aee358b4c6d0c8d6ccc36524c96e4c92ccf26e77e"
+dependencies = [
+ "num-integer",
+ "strength_reduce",
+]
+
+[[package]]
name = "unicode-ident"
version = "1.0.24"
source = "registry+https://github.com/rust-lang/crates.io-index"
diff --git a/Cargo.toml b/Cargo.toml
@@ -20,6 +20,7 @@ crate-type = ["cdylib", "rlib"]
# Uncomment the below line to disable the on-by-default VST3 feature to remove
# the GPL compatibility requirement
nih_plug = { git = "https://github.com/robbert-vdh/nih-plug.git", default-features = false, features = ["assert_process_allocs"] }
+rustfft = "6.4.1"
textplots = "0.8.7"
[profile.release]
diff --git a/README.md b/README.md
@@ -14,87 +14,107 @@ Purpose of this repository is creating an expressive wind synthesizer, like real
### Architecture
-#### Temporary oscillation
+#### Primary oscillation
This parts play a role of generating sounds like the reed on a saxophone or the lips on a trumpet.
<details>
<summary>TL;DR Derivation of the simulation formula</summary>
-## Derivation of the simulation formula
+```markdown
+## Physical Modeling and Discretization Process
-### 1. Continuous-Time Physical Equations
+This plugin simulates the sound generation mechanism of a reed instrument (or lip-reed instrument) by coupling a continuous-time mechanical oscillator with a non-linear fluid dynamics engine.
-#### Mechanical Oscillator (Reed Dynamics)
-
-The mechanical movement of the reed is modeled as a damped mass-spring system driven by an external fluid force $f(t)$:
-
-$$m \frac{d^2 x(t)}{dt^2} + r \frac{dx(t)}{dt} + k x(t) = f(t)$$
-
-Where:
+---
-* $m$: Effective mass of the reed.
-* $r$: Mechanical damping.
-* $k$: Stiffness (restoring force coefficient).
-* $x(t)$: Reed displacement ($x=0$ at rest, $x=2.0$ represents complete channel closure).
+### 1. Fluid Dynamics Discretization & Velocity Derivation
-#### Fluid Dynamics (Pressure-Velocity Relation)
+#### Continuous-Time Fluid Equation
-The airflow through the reed gap incorporates both the pressure drop across the orifice (Bernoulli's principle) and the acoustic/fluid inertia of the air mass within the channel:
+The pressure drop $P(t)$ across the orifice incorporates both the acoustic/fluid inertia of the air mass within the channel and Bernoulli's principle:
-$$\rho L \frac{dv_f(t)}{dt} + \frac{\rho}{4 g(t)^2} v_f(t)^2 = P(t)$$
+$$P(t) = \rho L \frac{dv_f(t)}{dt} + B(t) v_f(t)^2$$
Where:
* $\rho$: Air density ($1.2 \text{ kg/m}^3$).
* $L$: Effective length of the fluid column.
* $v_f(t)$: Fluid flow velocity.
-* $P(t)$: Driving breath pressure.
-* $g(t) = \max(2.0 - x(t), \epsilon)$: Dynamic aperture (gap width, clamped with a tiny $\epsilon$ to prevent division by zero).
+* $B(t) = \frac{\rho}{4 g(t)^2}$: Geometry-dependent flow resistance coefficient ($g(t)$ is the dynamic aperture).
+
+#### Bilinear Transform (Trapezoidal Integration)
+
+To discretize the derivative term, we apply the bilinear transform, which is mathematically equivalent to the trapezoidal rule. The derivative of fluid velocity at step $n$ is approximated as:
+
+$$\frac{dv_f(t)}{dt} \approx \frac{2}{T} (v_f[n] - v_f[n-1]) - \left.\frac{dv_f(t)}{dt}\right|_{n-1}$$
+
+Substituting the continuous fluid equation at step $n-1$ into the historic derivative term yields:
+
+$$\frac{dv_f(t)}{dt} \approx \frac{2}{T} (v_f[n] - v_f[n-1]) - \frac{1}{\rho L} \left( P[n-1] - B[n-1] v_f[n-1]^2 \right)$$
+
+Substituting this approximation back into the continuous-time equation at step $n$:
+
+$$P[n] = \rho L \left[ \frac{2}{T} (v_f[n] - v_f[n-1]) - \frac{1}{\rho L} \left( P[n-1] - B[n-1] v_f[n-1]^2 \right) \right] + B[n] v_f[n]^2$$
+
+Expanding and organizing the equation into a quadratic form with respect to the current velocity $v_f[n]$:
+
+$$B[n] v_f[n]^2 + \left( \frac{2\rho L}{T} \right) v_f[n] - \left[ P[n] + P[n-1] + \frac{2\rho L}{T} v_f[n-1] - B[n-1] v_f[n-1]^2 \right] = 0$$
+
+#### Analytical Solution for Discrete Fluid Velocity
+
+To solve for the physically valid (positive) root of this quadratic equation, let:
+
+* $A = \frac{2\rho L}{T}$
+* $C[n-1] = P[n] + P[n-1] + A v_f[n-1] - B[n-1] v_f[n-1]^2$
+
+Applying the quadratic formula explicitly determines the discrete fluid velocity $v_f[n]$ at the current time-step:
-#### Aeroelastic Coupling (Fluid Force)
+$$v_f[n] = \frac{-A + \sqrt{A^2 + 4 B[n] C[n-1]}}{2 B[n]}$$
-The aerodynamic force $f(t)$ acting on the surface of the reed depends on the dynamic pressure and the geometry of the channel:
+#### Mapping to Fluid Force $f[n]$
-$$f(t) = \pm \frac{1}{2} \rho v_f(t)^2 g(t)$$
+The calculated velocity $v_f[n]$ is immediately mapped to the aerodynamic force $f[n]$ acting on the reed surface based on the selected instrument mode:
-* **$\boldsymbol{+}$ (Positive Sign):** `SingleReed` Mode (Saxophone/Clarinet). The high velocity creates suction (Bernoulli effect) that pulls the reed toward closure.
-* **$\boldsymbol{-}$ (Negative Sign):** `LipReed` Mode (Trumpet/Brass). The pressure pushes the lips outward to open the channel.
+$$f[n] = \pm \frac{1}{2} \rho v_f[n]^2 g[n]$$
+
+* **$\boldsymbol{+}$ (Positive Sign):** `SingleReed` Mode (Saxophone). The high velocity creates suction (Bernoulli effect) that pulls the reed toward closure.
+* **$\boldsymbol{-}$ (Negative Sign):** `LipReed` Mode (Trumpet). The pressure pushes the lips outward to open the channel.
---
-### 2. Discretization via Bilinear Transform
+### 2. Mechanical Oscillator Discretization
-To compute this system inside a digital signal processor at sampling interval $T = 1 / f_s$, we map the continuous differential equations to discrete difference equations using the **Bilinear Transform (Tustin's Method)**.
+#### Continuous-Time Mechanical Equation
-The continuous derivative operator in the Laplace domain, $s$, is substituted by the discrete variable $z$ using the trapezoidal integration approximation:
+The movement of the mechanical reed is modeled as a damped mass-spring system driven by the derived fluid force $f(t)$:
-$$s \approx \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}$$
+$$m \frac{d^2 x(t)}{dt^2} + r \frac{dx(t)}{dt} + k x(t) = f(t)$$
-Applying this to the second-order derivative ($s^2$) yields:
+Where $m$ is the effective mass, $r$ is the mechanical damping, $k$ is the stiffness, and $x(t)$ is the displacement.
-$$s^2 \approx \frac{4}{T^2} \frac{1 - 2z^{-1} + z^{-2}}{1 + 2z^{-1} + z^{-2}}$$
+#### Bilinear Transform of the Oscillator
-#### Derivation of the Discrete Difference Equation
+We map the continuous differential system to the discrete $z$-domain by substituting the complex frequency $s$ via Tustin's method:
+
+$$s \approx \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}$$
-Substituting these into the mechanical transfer function $H(s) = \frac{X(s)}{F(s)} = \frac{1}{ms^2 + rs + k}$ gives:
+Applying this substitution to the second-order mechanical transfer function $H(s) = \frac{X(s)}{F(s)} = \frac{1}{ms^2 + rs + k}$ yields:
$$\frac{X(z)}{F(z)} = \frac{1}{m \left(\frac{4}{T^2} \frac{1 - 2z^{-1} + z^{-2}}{1 + 2z^{-1} + z^{-2}}\right) + r \left(\frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}\right) + k}$$
-Multiplying both the numerator and denominator by $(1 + 2z^{-1} + z^{-2})$ and grouping identical powers of $z^{-1}$, we clear the fraction fractions. To eliminate $1/T^2$ fractions and maximize numerical precision in single-precision floating-point math (`f32`), we multiply the entire equation by $T^2$:
+Multiplying both the numerator and denominator by $(1 + 2z^{-1} + z^{-2})$ and scaling the entire equation by $T^2$ to eliminate fractional sampling intervals ensures maximum numerical precision in single-precision floating-point math (`f32`):
$$\frac{X(z)}{F(z)} = \frac{T^2 (1 + 2z^{-1} + z^{-2})}{(4m + 2rT + kT^2) + (-8m + 2kT^2)z^{-1} + (4m - 2rT + kT^2)z^{-2}}$$
-Thus, we obtain the standard **Direct Form I** difference equation coefficients:
+This defines the standard **Direct Form I** difference equation coefficients:
-* $b_0 = T^2$
-* $b_1 = 2T^2$
-* $b_2 = T^2$
+* $b_0 = T^2, \quad b_1 = 2T^2, \quad b_2 = T^2$
* $a_0 = 4m + 2rT + kT^2$
* $a_1 = -8m + 2kT^2$
* $a_2 = 4m - 2rT + kT^2$
-The dynamic calculation of the exact displacement $x[n]$ at the current time-step is explicitly resolved as:
+The exact discrete displacement $x[n]$ at the current time-step is calculated as:
$$x[n] = \frac{b_0 f[n] + b_1 f[n-1] + b_2 f[n-2] - a_1 x[n-1] - a_2 x[n-2]}{a_0}$$
@@ -102,42 +122,131 @@ $$x[n] = \frac{b_0 f[n] + b_1 f[n-1] + b_2 f[n-2] - a_1 x[n-1] - a_2 x[n-2]}{a_0
### 3. Proof of Approximation Validity and Stability
-We prove that this discrete equation is a highly appropriate mathematical approximation of the continuous physical system $x(t)$ based on three criteria.
-
#### Proof A: Frequency Mapping Consistency
-The Bilinear Transform maps the entire continuous imaginary axis ($s = j\Omega$) onto the discrete unit circle ($z = e^{j\omega T}$). The mapping relationship is exactly:
+The Bilinear Transform maps the continuous imaginary axis ($s = j\Omega$) onto the discrete unit circle ($z = e^{j\omega T}$) via the exact relationship:
$$\Omega = \frac{2}{T} \tan\left(\frac{\omega T}{2}\right)$$
-For audio rates where the natural resonant frequency of the reed $\Omega_0 = \sqrt{k/m}$ satisfies $\Omega_0 \ll \frac{2}{T}$ (highly true since reed resonances are typically below $5\text{ kHz}$ and $T^{-1} = 44.1\text{ kHz}$), the Taylor expansion of the tangent function yields:
+For audio rates where the natural resonant frequency of the reed $\Omega_0 = \sqrt{k/m}$ satisfies $\Omega_0 \ll \frac{2}{T}$ (highly true since reed resonances are typically below $5\text{ kHz}$ and $T^{-1} = 44.1\text{ kHz}$), the Taylor expansion of the tangent function yields $\Omega \approx \omega$. This proves that the discrete resonance matches the continuous physical spectrum without severe high-frequency warping in the audible band.
-$$\Omega \approx \frac{2}{T} \left( \frac{\omega T}{2} \right) = \omega$$
+#### Proof B: Unconditional Numerical Stability
-This proves that in the audible band, the discrete frequency spectrum matches the physical continuous resonance behavior without severe warping.
+A physical reed system is passive and absorbs/dissipates energy via $r$. The continuous system poles lie in the Left-Half of the s-plane ($\text{Re}(s) = \sigma < 0$) because $m, r, k > 0$. Under the bilinear mapping:
-#### Proof B: Unconditional Numerical Stability (Passivity Preservation)
+$$|z|^2 = \left| \frac{1 + \frac{T}{2}s}{1 - \frac{T}{2}s} \right|^2 = \frac{(1 + \frac{T}{2}\sigma)^2 + (\frac{T}{2}\Omega)^2}{(1 - \frac{T}{2}\sigma)^2 + (\frac{T}{2}\Omega)^2}$$
-A physical reed is a passive system that absorbs and dissipates energy via $r$. For a system to be stable in the discrete domain, its poles must lie strictly inside the unit circle ($|z| < 1$).
+Since $\sigma < 0$, $(1 + \frac{T}{2}\sigma)^2 < (1 - \frac{T}{2}\sigma)^2$, mathematically guaranteeing $|z| < 1$.
-The continuous system poles lie in the Left-Half of the s-plane ($\text{Re}(s) < 0$) because $m, r, k > 0$. Under the bilinear mapping:
+> **Conclusion:** The system remains **unconditionally stable** regardless of real-time sampling rate modifications or aggressive parameter modulation ($m, k, r$ adjustments via bite intensity), eliminating numerical explosion risks common in forward-Euler methods.
+```
-$$z = \frac{1 + \frac{T}{2}s}{1 - \frac{T}{2}s}$$
+</details>
-Taking the magnitude squared when $\text{Re}(s) = \sigma < 0$:
+x(,f and v_f) formuler is
+
+$$x[n] = \frac{b_0 f[n] + b_1 f[n-1] + b_2 f[n-2] - a_1 x[n-1] - a_2 x[n-2]}{a_0}$$
-$$|z|^2 = \frac{(1 + \frac{T}{2}\sigma)^2 + (\frac{T}{2}\Omega)^2}{(1 - \frac{T}{2}\sigma)^2 + (\frac{T}{2}\Omega)^2}$$
+$$f[n] = \pm \frac{1}{2} \rho v_f[n]^2 g[n]$$
-Since $\sigma < 0$, it is mathematically guaranteed that $(1 + \frac{T}{2}\sigma)^2 < (1 - \frac{T}{2}\sigma)^2$, which proves $|z| < 1$.
+$$v_f[n] = \frac{-\alpha + \sqrt{\alpha^2 + 4 B[n] \Gamma[n-1]}}{2 B[n]}$$
-> **Conclusion:** The algorithm is **unconditionally stable** regardless of sample rate modifications or sudden physical parameter adjustments ($m, k, r$ modulation via bite intensity), eliminating numerical explosion risks inherent in explicit forward-Euler methods.
+#### Resonance Part
-#### Proof C: Order of Accuracy (Trapezoidal Match)
+<details>
+<summary>Physical Modeling of Boundary Dissipation in Open-Ended Acoustic Tubes</summary>
+
+This document provides the theoretical background, mathematical derivation, and discrete-time modeling of the acoustic reflection and energy dissipation at the open end (bell) of a resonant tube. Rather than using an empirical low-pass filter to attenuate high-frequency components, this model derives the high-frequency decay directly from the **physical mechanics of particle motion and nonlinear aerodynamic resistance** at the boundary.
+
+---
+
+## 1. Physical Phenomenon: Open-End Reflection and Turbulence
+
+When an acoustic wave propagating inside a tube reaches an open end, it transitions from a highly constrained, one-dimensional channel to a free, unconstrained three-dimensional space.
+
+1. **Inertial Overshoot:** As a high-pressure (condensation) wave arrives at the exit, the air particles are suddenly liberated and accelerate outward into the ambient atmosphere. Due to their mass (inertia), these particles overshoot, evacuation the region just inside the exit. This creates a local localized low-pressure (rarefaction) zone, which propagates back into the tube as a phase-inverted ($-1$) reflection.
+2. **Nonlinear Energy Dissipation ($v|v|$ Loss):** When air particles exhaust violently out of the tube orifice, the sharp edge causes flow separation, generating **local vortices and turbulence**. In fluid dynamics, this rapid change in cross-sectional area converts kinetic energy into unrecoverable thermal and acoustic radiation losses. The pressure drop $\Delta P$ across such an orifice is dominated by dynamic pressure, meaning it is proportional to the **square of the particle velocity** ($v^2$).
+
+---
+
+## 2. Continuous-Time Governing Equation
+
+We model the air plug at the open end as an effective acoustic mass (inertial slug) $M_e$ subject to a nonlinear aerodynamic resistance proportional to $v|v|$ to ensure that the damping force always opposes the direction of flow.
+
+Let $p_{total}(t)$ be the total acoustic pressure at the boundary, $S$ be the cross-sectional area of the tube, and $v(t)$ be the acoustic particle velocity. The force balance equation at the boundary $x = L$ is formulated as follows:
+
+$$S \cdot p_{total}(t) = M_e \frac{dv(t)}{dt} + r_{loss} v(t) |v(t)|$$
+
+Where:
+
+* $M_e = \rho_0 S \Delta L$ is the effective mass of the air plug (where $\Delta L \approx 0.613 \cdot \text{radius}$ is the open-end boundary correction length).
+* $r_{loss} = \frac{1}{2} \rho_0 S$ is the nonlinear loss coefficient derived from the dynamic pressure dissipation rate.
+* $\rho_0$ is the ambient air density.
+
+---
+
+## 3. Boundary Coupling with Digital Waveguide
+
+In a digital waveguide framework, the total pressure $p_{total}$ at the boundary is the sum of the incoming progressive wave (right-going wave $p^+$) and the reflected wave (left-going wave $p^-$):
+
+$$p_{total}(t) = p^+(t) + p^-(t)$$
+
+The particle velocity $v(t)$ is related to these traveling components via the characteristic acoustic impedance of the tube, $Z_0 = \frac{\rho_0 c}{S}$:
+
+$$v(t) = \frac{1}{Z_0} \left( p^+(t) - p^-(t) \right)$$
+
+By rewriting the reflected wave $p^-(t)$ in terms of the current velocity and the incident wave, we eliminate the algebraic loop:
+
+$$p^-(t) = p^+(t) - Z_0 v(t)$$
+
+Substituting this back into the total pressure equation yields:
+
+$$p_{total}(t) = 2p^+(t) - Z_0 v(t)$$
+
+---
+
+## 4. Discrete-Time Derivation via Backward Euler Discretization
+
+To solve the differential equation numerically without implicit algebraic loops, we apply a first-order Backward Euler finite difference approximation to the derivative term:
+
+$$\frac{dv(t)}{dt} \approx \frac{v[n] - v[n-1]}{T}$$
+
+Where $T = 1/f_s$ is the sampling period. Substituting the discrete derivative and the impedance-coupled total pressure into the continuous governing equation gives:
+
+$$S \left( 2p^+[n] - Z_0 v[n] \right) = M_e \left( \frac{v[n] - v[n-1]}{T} \right) + r_{loss} v[n] |v[n]|$$
+
+Assuming a positive outward flow ($v[n] > 0 \implies v[n]|v[n]| = v[n]^2$) during the primary exhaust cycle, we collect terms to form a standard quadratic equation with respect to the current velocity $v[n]$:
+
+$$r_{loss} v[n]^2 + \left( \frac{M_e}{T} + S Z_0 \right) v[n] - \left( 2S p^+[n] + \frac{M_e}{T} v[n-1] \right) = 0$$
+
+Let us define the discrete time-invariant and time-varying coefficients as:
+
+$$\begin{aligned}
+A_{bc} &= r_{loss} \\
+B_{bc} &= \frac{M_e}{T} + S Z_0 \\
+C_{bc}[n] &= 2S p^+[n] + \frac{M_e}{T} v[n-1]
+\end{aligned}$$
+
+The equation reduces to:
+
+$$A_{bc} v[n]^2 + B_{bc} v[n] - C_{bc}[n] = 0$$
+
+Applying the quadratic formula yields the explicit, stable algebraic solution for the particle velocity at the current sample step $n$:
+
+$$v[n] = \frac{-B_{bc} + \sqrt{B_{bc}^2 + 4 A_{bc} C_{bc}[n]}}{2 A_{bc}}$$
+
+Once $v[n]$ is evaluated, the true physical reflected wave $p^-[n]$ injected back into the delay line is determined deterministically:
+
+$$p^-[n] = p^+[n] - Z_0 v[n]$$
+
+---
-The bilinear transform is mathematically isomorphic to the trapezoidal integration rule. The local truncation error ($LTE$) of a trapezoidal approximation for a state vector $\mathbf{X}$ is bounded by:
+## 5. Why This Causes High-Frequency Decay
-$$LTE = \mathcal{O}(T^3)$$
+This mathematical formulation inherently explains why higher frequencies experience more severe damping without artificially implementing a digital filter:
-Accumulated over a global simulation window, the overall approximation error scales as $\mathcal{O}(T^2)$ (Second-order accurate). Compared to standard first-order Euler methods ($\mathcal{O}(T)$), this guarantees that high-frequency physical transcripts (such as rapid transients during attacks or lip-buzzing regimes) are preserved with minimal numerical artificial damping.
+1. **Velocity-Dependent Damping:** High-frequency components undergo rapid structural changes over time ($\frac{dv}{dt}$ is large), forcing higher transient particle velocities $v[n]$.
+2. **Quadratic Penalty:** Because the dissipative term scales with $v[n]^2$, these high-velocity, high-frequency transients encounter a nonlinearly magnified resistance compared to slow, low-frequency pressure oscillations.
+3. **Impedance Matching Shift:** At low velocities, the quadratic term vanishes, and the system matches an ideal open boundary, achieving an inversion coefficient close to $-1$. At high velocities, the effective boundary impedance changes due to $r_{loss}v^2$, allowing energy to escape into the environment as radiation or heat rather than reflecting back into the pipe.
</details>
diff --git a/src/lib.rs b/src/lib.rs
@@ -5,6 +5,190 @@ use std::sync::Arc;
// https://github.com/robbert-vdh/nih-plug/blob/master/plugins/examples/gain/src/lib.rs to get
// started
+const RHO: f32 = 1.2;
+const C_SPEED: f32 = 340.0;
+const RADIUS: f32 = 0.015;
+
+pub struct Bore {
+ pub delay_buffer: Vec<f32>,
+ pub pointer: usize,
+ pub z0: f32,
+ pub m_e: f32,
+ pub r_loss: f32,
+ pub v_prev: f32,
+}
+
+impl Bore {
+ pub fn new(target_freq: f32, sample_rate: f32, boundary_type: BoundaryType) -> Self {
+ let wavelength_factor = match boundary_type {
+ BoundaryType::NonlinearDissipative => 2.0,
+ BoundaryType::Fixed => 1.0,
+ };
+ let delay_samples = (sample_rate / (target_freq * wavelength_factor)).round() as usize;
+ let s = std::f32::consts::PI * RADIUS * RADIUS;
+ let z0 = (RHO * C_SPEED) / s;
+ let delta_l = 0.613 * RADIUS;
+ let m_e = RHO * s * delta_l;
+ let r_loss = 0.5 * RHO * s;
+
+ Self {
+ delay_buffer: vec![0.0; delay_samples],
+ pointer: 0,
+ z0,
+ m_e,
+ r_loss,
+ v_prev: 0.0,
+ }
+ }
+
+ pub fn set_frequency(
+ &mut self,
+ target_freq: f32,
+ sample_rate: f32,
+ boundary_type: BoundaryType,
+ ) {
+ let wavelength_factor = match boundary_type {
+ BoundaryType::NonlinearDissipative => 2.0,
+ BoundaryType::Fixed => 1.0,
+ };
+ let delay_samples = (sample_rate / (target_freq * wavelength_factor)).round() as usize;
+ if delay_samples != self.delay_buffer.len() && delay_samples > 0 {
+ self.delay_buffer.resize(delay_samples, 0.0);
+ if self.pointer >= delay_samples {
+ self.pointer = 0;
+ }
+ }
+ }
+
+ pub fn step(&mut self, p_plus_in: f32, t: f32, boundary_type: BoundaryType) -> f32 {
+ if self.delay_buffer.is_empty() {
+ return 0.0;
+ }
+
+ let p_minus_out = self.delay_buffer[self.pointer];
+
+ let p_minus_reflected = match boundary_type {
+ BoundaryType::NonlinearDissipative => {
+ let s = std::f32::consts::PI * RADIUS * RADIUS;
+ let b_bc = (self.m_e / t) + (s * self.z0);
+ let c_bc = 2.0 * s * p_plus_in + (self.m_e / t) * self.v_prev;
+
+ let discriminant = (b_bc * b_bc + 4.0 * self.r_loss * c_bc).max(0.0);
+ let v_current = (-b_bc + discriminant.sqrt()) / (2.0 * self.r_loss);
+
+ let res = p_plus_in - self.z0 * v_current;
+ self.v_prev = v_current;
+ res
+ }
+ BoundaryType::Fixed => p_plus_in,
+ };
+
+ self.delay_buffer[self.pointer] = p_minus_reflected;
+ self.pointer = (self.pointer + 1) % self.delay_buffer.len();
+
+ p_minus_out
+ }
+
+ pub fn reset(&mut self) {
+ for val in self.delay_buffer.iter_mut() {
+ *val = 0.0;
+ }
+ self.pointer = 0;
+ self.v_prev = 0.0;
+ }
+}
+
+pub struct Reed {
+ pub x_prev1: f32,
+ pub x_prev2: f32,
+ pub f_prev1: f32,
+ pub f_prev2: f32,
+ pub vf_prev1: f32,
+}
+
+impl Reed {
+ pub fn new() -> Self {
+ Self {
+ x_prev1: 0.0,
+ x_prev2: 0.0,
+ f_prev1: 0.0,
+ f_prev2: 0.0,
+ vf_prev1: 0.0,
+ }
+ }
+
+ pub fn step(
+ &mut self,
+ p_total: f32,
+ oscillation_type: OscillationType,
+ sample_rate: f32,
+ params: &AerothesisParams,
+ v_bite: f32,
+ ) -> (f32, f32) {
+ let t = 1.0 / sample_rate;
+ const EPS: f32 = 1e-5;
+
+ let gap_curr = (2.0 - self.x_prev1).max(EPS);
+ let b_curr = RHO / (4.0 * gap_curr * gap_curr);
+ let a_fluid = (RHO * params.reed_length.value()) / t;
+
+ let c_prev = p_total - b_curr * (self.vf_prev1 * self.vf_prev1);
+ let discriminant =
+ (a_fluid * a_fluid + 4.0 * b_curr * (a_fluid * self.vf_prev1 + c_prev)).max(0.0);
+ let vf_current = (-a_fluid + discriminant.sqrt()) / (2.0 * b_curr);
+
+ let mut f_current = if self.x_prev1 >= 2.0 {
+ 0.0
+ } else {
+ 0.5 * RHO * (vf_current * vf_current) * gap_curr
+ };
+
+ if oscillation_type == OscillationType::LipReed {
+ f_current = -f_current;
+ }
+
+ let m = params.base_mass.value() * (1.0 - params.bite_mass_scale.value() * v_bite);
+ let r = params.base_damping.value() * (1.0 + params.bite_damping_scale.value() * v_bite);
+ let k =
+ params.base_stiffness.value() * (1.0 + params.bite_stiffness_scale.value() * v_bite);
+
+ let b0 = t * t;
+ let b1 = 2.0 * t * t;
+ let b2 = t * t;
+
+ let a0 = 4.0 * m + 2.0 * r * t + k * t * t;
+ let a1 = -8.0 * m + 2.0 * k * t * t;
+ let a2 = 4.0 * m - 2.0 * r * t + k * t * t;
+
+ let mut x_n = (b0 * f_current + b1 * self.f_prev1 + b2 * self.f_prev2
+ - a1 * self.x_prev1
+ - a2 * self.x_prev2)
+ / a0;
+
+ if x_n >= 2.0 {
+ x_n = 2.0;
+ } else if x_n < 0.0 {
+ x_n = 0.0;
+ }
+
+ self.x_prev2 = self.x_prev1;
+ self.x_prev1 = x_n;
+ self.f_prev2 = self.f_prev1;
+ self.f_prev1 = f_current;
+ self.vf_prev1 = vf_current;
+
+ (x_n, vf_current)
+ }
+
+ pub fn reset(&mut self) {
+ self.x_prev1 = 0.0;
+ self.x_prev2 = 0.0;
+ self.f_prev1 = 0.0;
+ self.f_prev2 = 0.0;
+ self.vf_prev1 = 0.0;
+ }
+}
+
pub struct Aerothesis {
pub params: Arc<AerothesisParams>,
@@ -21,16 +205,29 @@ pub struct Aerothesis {
pub v_bite: f32,
pub v_fluid_prev: f32,
+
+ pub bore: Bore,
+ pub reed: Reed,
+ pub p_minus: f32,
+ pub current_frequency: f32,
}
#[derive(Enum, PartialEq, Clone, Copy)]
-pub enum InstrumentType {
+pub enum OscillationType {
#[name = "Single Reed"]
SingleReed,
- #[name = "Rip Reed"]
+ #[name = "Lip Reed"]
LipReed,
}
+#[derive(Enum, PartialEq, Clone, Copy)]
+pub enum BoundaryType {
+ #[name = "Nonlinear Dissipative"]
+ NonlinearDissipative,
+ #[name = "Fixed End"]
+ Fixed,
+}
+
#[derive(Params)]
pub struct AerothesisParams {
/// The parameter's ID is used to identify the parameter in the wrappred plugin API. As long as
@@ -40,8 +237,14 @@ pub struct AerothesisParams {
#[id = "gain"]
pub gain: FloatParam,
- #[id = "instrument_type"]
- pub instrument_type: EnumParam<InstrumentType>,
+ #[id = "output_gain"]
+ pub output_gain: FloatParam,
+
+ #[id = "oscillation_type"]
+ pub oscillation_type: EnumParam<OscillationType>,
+
+ #[id = "boundary_type"]
+ pub boundary_type: EnumParam<BoundaryType>,
#[id = "ReedLength"]
pub reed_length: FloatParam,
@@ -73,18 +276,26 @@ pub struct AerothesisParams {
impl Default for Aerothesis {
fn default() -> Self {
+ let sample_rate = 44100.0;
+ let target_freq = 220.0;
+ let params = Arc::new(AerothesisParams::default());
Self {
- params: Arc::new(AerothesisParams::default()),
+ bore: Bore::new(target_freq, sample_rate, params.boundary_type.value()),
+ params,
x_prev: 0.0,
x_prev2: 0.0,
v_prev: 0.0,
f_prev: 0.0,
f_prev2: 0.0,
- sample_rate: 44100.0,
+ sample_rate,
- v_breath: 0.1,
+ v_breath: 0.0,
v_bite: 0.0,
v_fluid_prev: 0.0,
+
+ reed: Reed::new(),
+ p_minus: 0.0,
+ current_frequency: target_freq,
}
}
}
@@ -116,7 +327,19 @@ impl Default for AerothesisParams {
.with_value_to_string(formatters::v2s_f32_gain_to_db(2))
.with_string_to_value(formatters::s2v_f32_gain_to_db()),
- instrument_type: EnumParam::new("Instrument Type", InstrumentType::SingleReed),
+ output_gain: FloatParam::new(
+ "Output Gain",
+ 0.01,
+ FloatRange::Skewed {
+ min: 0.0,
+ max: 1.0,
+ factor: 0.2,
+ },
+ ),
+
+ oscillation_type: EnumParam::new("Oscillation Type", OscillationType::SingleReed),
+
+ boundary_type: EnumParam::new("Boundary Type", BoundaryType::NonlinearDissipative),
reed_length: FloatParam::new(
"Reed Length",
@@ -356,12 +579,21 @@ impl Plugin for Aerothesis {
_context: &mut impl InitContext<Self>,
) -> bool {
self.sample_rate = buffer_config.sample_rate;
+ self.bore = Bore::new(
+ self.current_frequency,
+ self.sample_rate,
+ self.params.boundary_type.value(),
+ );
true
}
fn reset(&mut self) {
// Reset buffers and envelopes here. This can be called from the audio thread and may not
// allocate. You can remove this function if you do not need it.
+ self.bore.reset();
+ self.reed.reset();
+ self.p_minus = 0.0;
+ self.v_breath = 0.0;
}
fn process(
@@ -384,18 +616,52 @@ impl Plugin for Aerothesis {
self.v_bite = value;
}
}
+ NoteEvent::NoteOn { note, .. } => {
+ self.current_frequency = util::midi_note_to_freq(note);
+ self.bore.set_frequency(
+ self.current_frequency,
+ self.sample_rate,
+ self.params.boundary_type.value(),
+ );
+ }
_ => (),
}
}
+ let t = 1.0 / self.sample_rate;
+
for channel_samples in buffer.iter_samples() {
let gain = self.params.gain.smoothed.next();
- let x_n = self.step();
- self.v_fluid_prev = self.vf();
+ let p_total = self.v_breath + self.p_minus;
+ let (x_n, vf_n) = self.reed.step(
+ p_total,
+ self.params.oscillation_type.value(),
+ self.sample_rate,
+ &self.params,
+ self.v_bite,
+ );
+
+ let gap = (2.0 - x_n).max(1e-5);
+ let u = vf_n * gap;
+ let p_mouth = self.p_minus + self.bore.z0 * u;
+ let p_plus = p_mouth - self.p_minus;
+
+ self.p_minus = self.bore.step(p_plus, t, self.params.boundary_type.value());
+
+ let mut output = (p_mouth * self.params.output_gain.value()).clamp(-1.0, 1.0);
+
+ if self.v_breath < 1e-4 {
+ output = 0.0;
+ self.bore.reset();
+ self.reed.reset();
+ self.p_minus = 0.0;
+ }
+
+ let final_output = output * gain;
for sample in channel_samples {
- *sample = (x_n - self.equilibrium_offset()).clamp(-1.0, 1.0) * gain;
+ *sample = final_output;
}
}
@@ -414,6 +680,7 @@ impl ClapPlugin for Aerothesis {
ClapFeature::AudioEffect,
ClapFeature::Stereo,
ClapFeature::NoteDetector,
+ ClapFeature::Instrument,
];
}
diff --git a/src/main.rs b/src/main.rs
@@ -1,4 +1,5 @@
-use aerothesis::Aerothesis;
+use aerothesis::{Aerothesis, BoundaryType, OscillationType};
+use rustfft::{num_complex::Complex, FftPlanner};
use textplots::{Chart, Plot, Shape};
fn custom_step(plugin: &mut Aerothesis, m: f32, r: f32, k: f32) -> f32 {
@@ -59,35 +60,74 @@ fn custom_step(plugin: &mut Aerothesis, m: f32, r: f32, k: f32) -> f32 {
fn main() -> Result<(), Box<dyn std::error::Error>> {
let mut plugin = Aerothesis::default();
let sample_rate = 44100.0;
- let seconds = 0.02;
+ let seconds = 1.0;
let num_samples = (sample_rate * seconds) as usize;
let mut data = Vec::with_capacity(num_samples);
-
- // Balance for amplitude: m=0.0005, k=2000 (fn ~ 318Hz)
- let m = 0.0005;
- let r = 0.0001; // Extremely low damping
- let k = 2000.0;
+ let mut p_mouth_signal = Vec::with_capacity(num_samples);
plugin.sample_rate = sample_rate;
- println!(
- "Starting simulation for Large Amplitude (m={}, k={}, r={})...",
- m, k, r
- );
-
for i in 0..num_samples {
- // Very strong breath pressure
- plugin.v_breath = if i < 100 {
- (i as f32 / 100.0) * 100.0
+ plugin.v_breath = if i < 1000 {
+ (i as f32 / 1000.0) * 40.0
} else {
- 100.0
+ 40.0
};
- data.push((i as f32, custom_step(&mut plugin, m, r, k)));
+
+ let p_total = plugin.v_breath + plugin.p_minus;
+ let (x_n, vf_n) = plugin.reed.step(
+ p_total,
+ plugin.params.oscillation_type.value(),
+ plugin.sample_rate,
+ &plugin.params,
+ plugin.v_bite,
+ );
+
+ let gap = (2.0 - x_n).max(1e-5);
+ let u = vf_n * gap;
+ let p_mouth = plugin.p_minus + plugin.bore.z0 * u;
+ let p_plus = p_mouth - plugin.p_minus;
+
+ plugin.p_minus = plugin.bore.step(
+ p_plus,
+ 1.0 / sample_rate,
+ BoundaryType::NonlinearDissipative,
+ );
+
+ if i < (sample_rate * 0.05) as usize {
+ data.push((i as f32, p_mouth));
+ }
+ p_mouth_signal.push(p_mouth);
}
- Chart::new(180, 60, 0.0, num_samples as f32)
+ Chart::new(180, 60, 0.0, data.len() as f32)
.lineplot(&Shape::Lines(&data))
.display();
+ let fft_len = p_mouth_signal.len().next_power_of_two();
+ let mut planner = FftPlanner::new();
+ let fft = planner.plan_fft_forward(fft_len);
+
+ let mut buffer: Vec<Complex<f32>> = p_mouth_signal
+ .iter()
+ .map(|&p| Complex { re: p, im: 0.0 })
+ .collect();
+ buffer.resize(fft_len, Complex { re: 0.0, im: 0.0 });
+
+ fft.process(&mut buffer);
+
+ let mut spectrum_data = Vec::with_capacity(fft_len / 2);
+ for (i, complex) in buffer.iter().enumerate().take(fft_len / 2) {
+ let freq = (i as f32 * sample_rate) / fft_len as f32;
+ if freq > 1500.0 {
+ break;
+ }
+ spectrum_data.push((freq, complex.norm()));
+ }
+
+ Chart::new(180, 60, 0.0, 1500.0)
+ .lineplot(&Shape::Lines(&spectrum_data))
+ .display();
+
Ok(())
}
