commit 09a4efe8947563cb18b2cdf659abceab8ea62a99
parent 6bdb5e079c6e3fabb5f667e42f4eced3bfc4d7b8
Author: minerva-jupiter <ryouturn@gmail.com>
Date: Fri, 19 Jun 2026 12:15:28 +0900
feat(engine): refactor acoustic modeling to use discrete delay resonance
- Remove Bore and Reed structs in favor of a centralized resonance implementation using x_history and VecDeque.
- Update Aerothesis to support OpenPipe and ClosedPipe resonance types with configurable decay.
- Cleanup documentation and simulation code in main.rs.
Diffstat:
| M | README.md | | | 94 | ------------------------------------------------------------------------------- |
| M | src/lib.rs | | | 357 | +++++++++++++++++++------------------------------------------------------------ |
| M | src/main.rs | | | 117 | +++++++++++++++++++------------------------------------------------------------- |
3 files changed, 110 insertions(+), 458 deletions(-)
diff --git a/README.md b/README.md
@@ -150,97 +150,3 @@ $$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]}$$
#### Resonance Part
-
-<details>
-<summary>Physical Modeling of Boundary Dissipation in Open-Ended Acoustic Tubes</summary>
-
-## 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]$$
-
----
-
-## 5. Why This Causes High-Frequency Decay
-
-This mathematical formulation inherently explains why higher frequencies experience more severe damping without artificially implementing a digital filter:
-
-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
@@ -1,194 +1,10 @@
use nih_plug::prelude::*;
-use std::sync::Arc;
+use std::{collections::VecDeque, sync::Arc};
// This is a shortened version of the gain example with most comments removed, check out
// 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>,
@@ -206,26 +22,25 @@ pub struct Aerothesis {
pub v_fluid_prev: f32,
- pub bore: Bore,
- pub reed: Reed,
- pub p_minus: f32,
- pub current_frequency: f32,
+ pub x_history: VecDeque<f32>,
+
+ pub note_frequency: f32,
}
#[derive(Enum, PartialEq, Clone, Copy)]
-pub enum OscillationType {
+pub enum InstrumentType {
#[name = "Single Reed"]
SingleReed,
- #[name = "Lip Reed"]
+ #[name = "Rip Reed"]
LipReed,
}
#[derive(Enum, PartialEq, Clone, Copy)]
-pub enum BoundaryType {
- #[name = "Nonlinear Dissipative"]
- NonlinearDissipative,
- #[name = "Fixed End"]
- Fixed,
+pub enum ResonanceType {
+ #[name = "Open Pipe"]
+ OpenPipe,
+ #[name = "Closed Pipe"]
+ ClosedPipe,
}
#[derive(Params)]
@@ -237,14 +52,8 @@ pub struct AerothesisParams {
#[id = "gain"]
pub gain: FloatParam,
- #[id = "output_gain"]
- pub output_gain: FloatParam,
-
- #[id = "oscillation_type"]
- pub oscillation_type: EnumParam<OscillationType>,
-
- #[id = "boundary_type"]
- pub boundary_type: EnumParam<BoundaryType>,
+ #[id = "instrument_type"]
+ pub instrument_type: EnumParam<InstrumentType>,
#[id = "ReedLength"]
pub reed_length: FloatParam,
@@ -272,30 +81,32 @@ pub struct AerothesisParams {
pub breath_cc: IntParam,
#[id = "bite_cc"]
pub bite_cc: IntParam,
+
+ #[id = "resonance_type"]
+ pub resonance_type: EnumParam<ResonanceType>,
+
+ #[id = "resonance_decay"]
+ pub resonance_decay: FloatParam,
}
impl Default for Aerothesis {
fn default() -> Self {
- let sample_rate = 44100.0;
- let target_freq = 220.0;
- let params = Arc::new(AerothesisParams::default());
Self {
- bore: Bore::new(target_freq, sample_rate, params.boundary_type.value()),
- params,
+ params: Arc::new(AerothesisParams::default()),
x_prev: 0.0,
x_prev2: 0.0,
v_prev: 0.0,
f_prev: 0.0,
f_prev2: 0.0,
- sample_rate,
+ sample_rate: 44100.0,
- v_breath: 0.0,
+ v_breath: 0.1,
v_bite: 0.0,
v_fluid_prev: 0.0,
- reed: Reed::new(),
- p_minus: 0.0,
- current_frequency: target_freq,
+ x_history: VecDeque::new(),
+
+ note_frequency: 0.0,
}
}
}
@@ -327,19 +138,7 @@ 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()),
- 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),
+ instrument_type: EnumParam::new("Instrument Type", InstrumentType::SingleReed),
reed_length: FloatParam::new(
"Reed Length",
@@ -420,6 +219,17 @@ impl Default for AerothesisParams {
breath_cc: IntParam::new("Breath CC", 2, IntRange::Linear { min: 0, max: 127 }),
bite_cc: IntParam::new("Bite CC", 11, IntRange::Linear { min: 0, max: 127 }),
+
+ resonance_type: EnumParam::new("Resonance Type", ResonanceType::OpenPipe),
+ resonance_decay: FloatParam::new(
+ "Resonance Decay",
+ 0.9,
+ FloatRange::Skewed {
+ min: 0.0,
+ max: 1.0,
+ factor: 0.8,
+ },
+ ),
}
}
}
@@ -465,6 +275,7 @@ impl Aerothesis {
let b_prev = RHO / (4.0 * (gap_prev * gap_prev));
let c_prev = self.v_breath - b_prev * (self.v_fluid_prev * self.v_fluid_prev);
+ // Current gap is also based on x_prev in this discrete model for stability
let gap_curr = (2.0 - self.x_prev).clamp(EPS, 2.0);
let b_curr = RHO / (4.0 * (gap_curr * gap_curr));
@@ -519,6 +330,29 @@ impl Aerothesis {
let t = 1.0 / self.sample_rate;
(2.0 / t) * (x - self.x_prev) - self.v_prev
}
+
+ 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.x_history.len() as f32 {
+ 0.0
+ } else {
+ let decay: f32 = if self.params.resonance_type.value() == ResonanceType::OpenPipe {
+ 1.0
+ } else {
+ -1.0
+ } * self.params.resonance_decay.value();
+ let x_delay = self.x_history.pop_front().unwrap_or(0.0);
+ decay * x_delay
+ };
+
+ let x_current = x_oscillator + resonance;
+ self.x_history.push_back(x_current);
+
+ x_current
+ }
+
fn equilibrium_offset(&self) -> f32 {
let f = self.f();
let k = self.k();
@@ -528,6 +362,16 @@ impl Aerothesis {
0.0
}
}
+ 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.x_history.iter().sum::<f32>() / self.x_history.len() as f32
+ }
}
impl Plugin for Aerothesis {
@@ -578,21 +422,13 @@ 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;
+ self.x_history.clear();
}
fn process(
@@ -615,52 +451,26 @@ 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(),
- );
+ NoteEvent::NoteOn {
+ timing: _,
+ voice_id: _,
+ channel: _,
+ note,
+ velocity: _,
+ } => {
+ self.reset();
+ self.note_frequency = util::midi_note_to_freq(note);
}
_ => (),
}
}
- let t = 1.0 / self.sample_rate;
-
for channel_samples in buffer.iter_samples() {
let gain = self.params.gain.smoothed.next();
-
- 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;
+ let x_current = self.resonance() - self.avg_x_history();
for sample in channel_samples {
- *sample = final_output;
+ *sample = (x_current * gain).clamp(-1.0, 1.0);
}
}
@@ -679,7 +489,6 @@ impl ClapPlugin for Aerothesis {
ClapFeature::AudioEffect,
ClapFeature::Stereo,
ClapFeature::NoteDetector,
- ClapFeature::Instrument,
];
}
diff --git a/src/main.rs b/src/main.rs
@@ -1,117 +1,53 @@
-use aerothesis::{Aerothesis, BoundaryType, OscillationType};
+use aerothesis::Aerothesis;
+use nih_plug::prelude::util;
use rustfft::{num_complex::Complex, FftPlanner};
use textplots::{Chart, Plot, Shape};
-fn custom_step(plugin: &mut Aerothesis, m: f32, r: f32, k: f32) -> f32 {
- let t = 1.0 / plugin.sample_rate;
- const EPS: f32 = 1e-5;
- const RHO: f32 = 1.2;
-
- // vf logic (Current flow velocity)
- let vf = {
- let a_fluid = (RHO * plugin.params.reed_length.value()) / t;
- let gap_prev = (2.0 - plugin.x_prev).clamp(EPS, 2.0);
- let b_prev = RHO / (4.0 * (gap_prev * gap_prev));
- let c_prev = plugin.v_breath - b_prev * (plugin.v_fluid_prev * plugin.v_fluid_prev);
-
- let gap_curr = (2.0 - plugin.x_prev).clamp(EPS, 2.0);
- let b_curr = RHO / (4.0 * (gap_curr * gap_curr));
-
- if gap_curr <= EPS {
- 0.0
- } else {
- let disc = (a_fluid * a_fluid
- + 4.0 * b_curr * (a_fluid * plugin.v_fluid_prev + c_prev))
- .max(0.0);
- (-a_fluid + disc.sqrt()) / (2.0 * b_curr)
- }
- };
-
- // f logic (Force)
- let f = {
- let gap_curr = (2.0 - plugin.x_prev).clamp(EPS, 2.0);
- 0.5 * RHO * (vf * vf) * gap_curr
- };
-
- // x logic (Displacement)
- let x_n = {
- 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;
- ((b0 * f + b1 * plugin.f_prev + b2 * plugin.f_prev2
- - a1 * plugin.x_prev
- - a2 * plugin.x_prev2)
- / a0)
- .clamp(0.0, 2.0)
- };
-
- plugin.x_prev2 = plugin.x_prev;
- plugin.x_prev = x_n;
- plugin.f_prev2 = plugin.f_prev;
- plugin.f_prev = f;
- plugin.v_fluid_prev = vf;
-
- x_n
-}
-
fn main() -> Result<(), Box<dyn std::error::Error>> {
let mut plugin = Aerothesis::default();
let sample_rate = 44100.0;
- let seconds = 1.0;
+ let seconds = 0.5;
let num_samples = (sample_rate * seconds) as usize;
- let mut data = Vec::with_capacity(num_samples);
- let mut p_mouth_signal = Vec::with_capacity(num_samples);
plugin.sample_rate = sample_rate;
+ plugin.note_frequency = util::midi_note_to_freq(48); // C3 (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);
for i in 0..num_samples {
- plugin.v_breath = if i < 1000 {
- (i as f32 / 1000.0) * 40.0
+ // Simple attack envelope for breath pressure
+ plugin.v_breath = if i < 2000 {
+ (i as f32 / 2000.0) * 100.0
} else {
- 40.0
+ 100.0
};
- 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,
- );
+ let sample = plugin.resonance();
+ // Collect first 50ms for waveform plot
if i < (sample_rate * 0.05) as usize {
- data.push((i as f32, p_mouth));
+ data.push((i as f32, sample));
}
- p_mouth_signal.push(p_mouth);
+ signal.push(sample);
}
+ println!("--- Waveform (first 50ms) ---");
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 fft_len = 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();
+ let mut buffer: Vec<Complex<f32>> =
+ signal.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);
@@ -119,13 +55,14 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
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 {
+ if freq > 2000.0 {
break;
}
spectrum_data.push((freq, complex.norm()));
}
- Chart::new(180, 60, 0.0, 1500.0)
+ println!("\n--- Spectrum (0 - 2000Hz) ---");
+ Chart::new(180, 60, 0.0, 2000.0)
.lineplot(&Shape::Lines(&spectrum_data))
.display();
