commit 9b4a1596b9f00ee89a2daf7cfcba19015e1683e4
parent a5c7adfb95954021d459f80873c69ce786f8c74a
Author: minerva-jupiter <ryouturn@gmail.com>
Date: Sat, 20 Jun 2026 22:00:51 +0900
feat: implement sign-preserving cubic damping and refine resonance model
- Update physical model documentation to reflect the new sign-preserving cubic damping implementation.
- Refactor `resonance` and `displacement` logic in `Aerothesis` to prevent DC offset and signal rectification.
- Adjust resonance delay calculations for open and closed pipes.
- Add sample-rate-based buffer clearing to improve stability when no note is active.
- Update simulation parameters in `main.rs` and improve the visualization resolution.
Diffstat:
| M | README.md | | | 27 | ++++++++++++++++----------- |
| M | src/lib.rs | | | 67 | ++++++++++++++++++++++++++++++++++++++++++------------------------- |
| M | src/main.rs | | | 43 | +++++++++++++++++++++++++++---------------- |
3 files changed, 85 insertions(+), 52 deletions(-)
diff --git a/README.md b/README.md
@@ -10,13 +10,13 @@ cargo xtask bundle aerothesis --release
## Design
-Purpose of this repository is creating an expressive wind synthesizer, like real trumpets, saxophones and other instruments.
+The purpose of this repository is to create an expressive wind synthesizer, simulating real trumpets, saxophones, and other wind instruments.
### Architecture
#### Primary oscillation
-This parts play a role of generating sounds like the reed on a saxophone or the lips on a trumpet.
+This part plays the role of generating sound, simulating the reed on a saxophone or the lips on a trumpet.
<details>
<summary>TL;DR Derivation of the simulation formula</summary>
@@ -141,13 +141,13 @@ Since $\sigma < 0$, $(1 + \frac{T}{2}\sigma)^2 < (1 - \frac{T}{2}\sigma)^2$, mat
</details>
-x(,f and v_f) formuler is
+The formula for $x[n]$, $f[n]$, and $v_f[n]$ 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}$$
$$f[n] = \pm \frac{1}{2} \rho v_f[n]^2 g[n]$$
-$$v_f[n] = \frac{-\alpha + \sqrt{\alpha^2 + 4 B[n] \Gamma[n-1]}}{2 B[n]}$$
+$$v_f[n] = \frac{-A + \sqrt{A^2 + 4 B[n] C[n-1]}}{2 B[n]}$$
#### Resonance Part
@@ -158,9 +158,9 @@ Rather than simulating wave reflection through complex fluid dynamics (changes i
#### 1. Damping Mechanism
-Energy in an acoustic system is proportional to the square of the time derivative of displacement ($(\partial x / \partial t)^2$). We apply a damping constant $a$ to this derivative. This effectively attenuates higher-frequency components, as their energy dissipates faster than lower-frequency components. Given an input displacement $x_n$, a delayed resonant displacement $x_{\text{resonance}}$, and the total previous displacement $x_{\text{prev}}$, the system state is updated as:
+Energy in an acoustic system is proportional to the square of the time derivative of displacement ($(\partial x / \partial t)^2$). We apply a damping constant $a$ to this derivative. To preserve the sign of the wave (preventing signal rectification and DC offset), the damping is implemented as a sign-preserving cubic non-linearity. Given an input displacement $x_{\text{in}}[n]$, a delayed resonant displacement $x_{\text{resonance}}[n]$, and the total previous displacement $x[n-1]$, the system state $x[n]$ is updated as:
-$$x = a (x_{\text{prev}} - (x_n + x_{\text{resonance}}))^2$$
+$$x[n] = (x_{\text{in}}[n] + x_{\text{resonance}}[n]) \cdot a \cdot (x[n-1] - (x_{\text{in}}[n] + x_{\text{resonance}}[n]))^2$$
#### 2. Physical Validity (D’Alembert’s Solution)
@@ -184,8 +184,8 @@ Consequently, calculating the resonance by multiplying the previously delayed di
Given a note frequency $f$ and the speed of sound $c$, the wavelength $\lambda$ is defined as $\lambda = c/f$.
-* **Open Pipe:** $\lambda = 2L \implies \text{round-trip time} = 2L/c = 1/f$.
-* **Closed Pipe:** $\lambda = 4L \implies \text{round-trip time} = 4L/c = 2/f$.
+* **Open Pipe:** $\lambda = 2L \implies \text{round-trip time} = 2L/c = \frac{2}{c} \frac{c}{2f} = \frac{1}{f}$.
+* **Closed Pipe:** $\lambda = 4L \implies \text{round-trip time} = 2L/c = \frac{2}{c} \frac{c}{4f} = \frac{1}{2f}$.
Thus, the required delay samples can be derived directly from the frequency $f$ and sample rate $fs$ without needing explicit values for tube length $L$ or sound speed $c$.
@@ -197,18 +197,23 @@ The core simulation is based on a displacement-driven delay-line model, where th
1. System Update Equation
-The total displacement $x[n]$ is calculated as a damped non-linear function of the input and the delayed resonant state. Given a damping constant $a$ ($0 < a \le 1$):
+The total displacement $x[n]$ is calculated as a damped non-linear function of the input and the delayed resonant state. To preserve the sign of the displacement wave and avoid DC rectification, a sign-preserving cubic function is used. Given a damping constant $a$ ($0 < a \le 1$):
-$$x[n] = a \cdot \left( x[n-1] - (x_{\text{in}}[n] + x_{\text{resonance}}[n]) \right)^2$$
+$$x[n] = (x_{\text{in}}[n] + x_{\text{resonance}}[n]) \cdot a \cdot \left( x[n-1] - (x_{\text{in}}[n] + x_{\text{resonance}}[n]) \right)^2$$
Where $x[n-1]$ represents the previous total displacement, capturing the system's memory.
2. Resonant Feedback (Delay and Reflection)
-The resonant component $x_{\text{resonance}}$ is the delayed state derived from the pipe's boundary conditions. Given a delay buffer $D$ of length $T$ (where $T = f_s / f$), the resonance is defined by the reflection coefficient $R$:
+The resonant component $x_{\text{resonance}}$ is the delayed state derived from the pipe's boundary conditions. Given a delay buffer $D$ of length $T$, the resonance is defined by the reflection coefficient $R$:
$$x_{\text{resonance}}[n] = R \cdot \text{buffer}[n - T]$$
+Where the round-trip delay length $T$ in samples is defined as:
+* **Open Pipes:** $T = \frac{f_s}{f}$
+* **Closed Pipes:** $T = \frac{f_s}{2f}$
+
+
* **For Open Pipes (Open-Open):**
* Reflection occurs twice per round-trip with a phase inversion, resulting in $R = 1$ (net phase preserved).
diff --git a/src/lib.rs b/src/lib.rs
@@ -22,9 +22,11 @@ pub struct Aerothesis {
pub v_fluid_prev: f32,
- pub displacement_history: VecDeque<f32>,
- pub displacement_prev: f32,
+ pub x_history: VecDeque<f32>,
+
pub note_frequency: f32,
+
+ pub displacement_prev: f32,
}
#[derive(Enum, PartialEq, Clone, Copy)]
@@ -104,9 +106,11 @@ impl Default for Aerothesis {
v_bite: 0.0,
v_fluid_prev: 0.0,
- displacement_history: VecDeque::new(),
- displacement_prev: 0.0,
+ x_history: VecDeque::new(),
+
note_frequency: 0.0,
+
+ displacement_prev: 0.0,
}
}
}
@@ -223,7 +227,7 @@ impl Default for AerothesisParams {
resonance_type: EnumParam::new("Resonance Type", ResonanceType::OpenPipe),
resonance_decay: FloatParam::new(
"Resonance Decay",
- 0.01,
+ 0.9,
FloatRange::Skewed {
min: 0.0,
max: 1.0,
@@ -332,31 +336,40 @@ impl Aerothesis {
}
pub fn resonance(&mut self) -> f32 {
- let x_n = self.step();
- let x_oscillator = x_n - self.equilibrium_offset();
-
- let resonance = if self.resonance_delay_samples() > self.displacement_history.len() as f32 {
+ if self.resonance_delay_samples() > self.x_history.len() as f32 {
0.0
} else {
+ if self.resonance_delay_samples() < self.x_history.len() as f32 {
+ self.x_history
+ .truncate(self.resonance_delay_samples() as usize);
+ }
let decay: f32 = if self.params.resonance_type.value() == ResonanceType::OpenPipe {
1.0
} else {
-1.0
};
- let x_delay = self.displacement_history.pop_front().unwrap_or(0.0);
+ let x_delay = self.x_history.pop_back().unwrap_or(0.0);
decay * x_delay
- };
+ }
+ }
- let x_nondamping = x_oscillator + resonance;
+ pub fn displacement(&mut self) -> f32 {
+ let x_n = self.step();
+ let x_oscillator = x_n - self.equilibrium_offset();
+
+ let resonance = self.resonance();
- let x_current = x_nondamping
- * (1.0 - self.params.resonance_decay.value())
- * (self.displacement_prev - x_nondamping)
- * (self.displacement_prev - x_nondamping);
+ let x_current = x_oscillator + resonance;
- self.displacement_history.push_back(x_current);
+ let displacement = x_current
+ * (self.params.resonance_decay.value()
+ * (x_current - self.displacement_prev)
+ * (x_current - self.displacement_prev))
+ .clamp(0.0, 1.0);
+ self.displacement_prev = displacement;
- x_current
+ self.x_history.push_front(displacement);
+ x_oscillator
}
fn equilibrium_offset(&self) -> f32 {
@@ -368,15 +381,15 @@ impl Aerothesis {
0.0
}
}
- fn resonance_delay_samples(&self) -> f32 {
+ pub fn resonance_delay_samples(&self) -> f32 {
if self.params.resonance_type.value() == ResonanceType::OpenPipe {
self.sample_rate / self.note_frequency
} else {
self.sample_rate / 2.0 / self.note_frequency
}
}
- fn avg_x_history(&self) -> f32 {
- self.displacement_history.iter().sum::<f32>() / self.displacement_history.len() as f32
+ pub fn avg_x_history(&self) -> f32 {
+ self.x_history.iter().sum::<f32>() / self.x_history.len() as f32
}
}
@@ -434,8 +447,7 @@ impl Plugin for Aerothesis {
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.displacement_history.clear();
+ self.x_history.clear();
}
fn process(
@@ -474,10 +486,15 @@ impl Plugin for Aerothesis {
for channel_samples in buffer.iter_samples() {
let gain = self.params.gain.smoothed.next();
- let x_current = self.resonance() - self.avg_x_history();
+ let x_current = self.displacement() - self.avg_x_history();
for sample in channel_samples {
- *sample = (x_current * gain).clamp(-1.0, 1.0);
+ if self.note_frequency == 0.0 {
+ self.x_history.clear();
+ *sample = 0.0;
+ } else {
+ *sample = (x_current * gain).clamp(-1.0, 1.0);
+ }
}
}
diff --git a/src/main.rs b/src/main.rs
@@ -6,21 +6,28 @@ use textplots::{Chart, Plot, Shape};
fn main() -> Result<(), Box<dyn std::error::Error>> {
let mut plugin = Aerothesis::default();
let sample_rate = 44100.0;
- let seconds = 0.5;
+ let seconds = 1.0;
let num_samples = (sample_rate * seconds) as usize;
plugin.sample_rate = sample_rate;
- plugin.note_frequency = util::midi_note_to_freq(48); // C3 (u8)
+ // plugin.note_frequency = util::midi_note_to_freq(48); // C3 (u8)
+ plugin.note_frequency = util::midi_note_to_freq(54); // C4 (u8)
// For simulation in main.rs, we use the default parameters from AerothesisParams::default()
// because nih-plug parameters are designed to be managed by a host and don't have
// simple setter methods for plain values without a ParamSetter context.
// Default resonance: OpenPipe, Decay: 0.9
- let mut data = Vec::with_capacity(num_samples);
- let mut signal = Vec::with_capacity(num_samples);
+ let mut displacements = Vec::with_capacity(num_samples);
+ let mut resonances = Vec::with_capacity(num_samples);
for i in 0..num_samples {
+ // let x_current = self.resonance() - self.avg_x_history();
+
+ // for sample in channel_samples {
+ // *sample = (x_current * gain).clamp(-1.0, 1.0);
+ // }
+
// Simple attack envelope for breath pressure
plugin.v_breath = if i < 2000 {
(i as f32 / 2000.0) * 100.0
@@ -28,26 +35,30 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
100.0
};
- let sample = plugin.resonance();
-
- // Collect first 50ms for waveform plot
- if i < (sample_rate * 0.05) as usize {
- data.push((i as f32, sample));
- }
- signal.push(sample);
+ let sample = plugin.displacement() - plugin.avg_x_history();
+ displacements.push(sample);
+ resonances.push(plugin.resonance());
}
- println!("--- Waveform (first 50ms) ---");
- Chart::new(180, 60, 0.0, data.len() as f32)
+ let len = (sample_rate * 0.05) as usize;
+
+ let data: Vec<(f32, f32)> = (0..len)
+ .map(|i| (i as f32, displacements[i + (sample_rate * 0.1) as usize]))
+ .collect();
+
+ println!("--- Waveform ---");
+ Chart::new(360, 60, 0.0, len as f32)
.lineplot(&Shape::Lines(&data))
.display();
- let fft_len = signal.len().next_power_of_two();
+ let fft_len = displacements.len().next_power_of_two();
let mut planner = FftPlanner::new();
let fft = planner.plan_fft_forward(fft_len);
- let mut buffer: Vec<Complex<f32>> =
- signal.iter().map(|&s| Complex { re: s, im: 0.0 }).collect();
+ let mut buffer: Vec<Complex<f32>> = displacements
+ .iter()
+ .map(|&s| Complex { re: s, im: 0.0 })
+ .collect();
buffer.resize(fft_len, Complex { re: 0.0, im: 0.0 });
fft.process(&mut buffer);
