diff --git a/CMakeLists.txt b/CMakeLists.txt index 17df385..f9a04fe 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -56,6 +56,7 @@ add_library(sonobulus_core STATIC src/synth/Synth.cpp src/synth/Voice.cpp src/synth/Instrument.cpp + src/synth/Instruments/PianoString.cpp ) target_link_libraries(sonobulus_core PRIVATE Qt6::Core diff --git a/scripts/string_model_fdm.py b/scripts/string_model_fdm.py index 5ee84eb..868c144 100644 --- a/scripts/string_model_fdm.py +++ b/scripts/string_model_fdm.py @@ -3,9 +3,10 @@ import math import numpy as np import matplotlib.pyplot as plt import sounddevice as sd +import time sample_rate = 44100 -seconds = 5 +seconds = 10 N = 100 # number of string segments I = int(sample_rate * seconds) # number of samples to simulate @@ -19,23 +20,29 @@ T = 1200 # string tension, N c = math.sqrt(T/mu) # transverse wave velocity kappa = 0.001 # stiffness coefficient sigma = 0.5 # damping coefficient -L = 0.5 # length of string +L = 1.0 # length of string strike_position = 0.2 # x of impulse location impulse_width = 0.02 # x of impulse width impulse_velocity = 1000.0 # x/t of impulse magnitude sample_position = 0.1 # percentage along L of sampling for audio +f_0 = c / (2*L) # fundamental frequency of a non-stiff string +f_1 = f_0 * math.sqrt(1 + kappa) # fundamental frequency of the stiff string +print("fundamental frequency =", f_1) + dx = L / N # delta x dt = 1/sample_rate #dt = 0.2 * dx / c if(dx**2 < (c*dt)**2 + 4*(kappa*dt/(dx**2))**2): - print("warning: possibly unstable due to not enough segments, increase N") + print("warning: possibly unstable, increase sample rate or decrease string segments") # derived constants r1 = c * dt/dx r2 = (c * dt/dx) ** 2 s1 = kappa * dt/dx**2 s2 = (kappa * dt/dx**2) ** 2 +a1 = 2 - 2*sigma*dt +a2 = 2*sigma*dt - 1 # string grid x = np.linspace(0, L, N + 1) # linspace my beloved @@ -59,16 +66,16 @@ def show_plot(): plt.grid() plt.show() +start_time = time.perf_counter() + +def applyImpulse(): + for n in range(2, N-2): + y_current[n] = y_current[n] + dt*v0[n] + # first iteration for n in range(2, N-2): - y_current[n] = dt * v0[n] + 0.5 * r1**2 * (y_last[n-1] - 2*y_last[n] + y_last[n+1]) - # y_xx = y_last[n+1] - 2*y_last[n] + y_last[n-1] - # y_xxxx = y_last[n-2] - 4*y_last[n-1] + 6*y_last[n] - 4*y_last[n+1] + y_last[n+2] - - # term1 = (2 - 2*sigma*dt) * y_current[n] - # term2 = (2*sigma*dt - 1) * y_last[n] - # term3 = r2 * y_xx - s2 * y_xxxx - # y_current[n] = dt * v0[n] + 0.5 * r1**2 *(term1 + term2 + term3) + y_current[n] = 0.5 * r1**2 * (y_last[n-1] - 2*y_last[n] + y_last[n+1]) +applyImpulse() n_sample = int(sample_position * L * N) y_sample[0] = y_last[n_sample] @@ -78,17 +85,11 @@ y_sample[1] = y_current[n_sample] for i in range(2, I): for n in range(2, N-2): - # simple wave equation - # y_next[n] = 2*y_current[n] - y_last[n] + r**2 * (y_current[n-1] - 2*y_current[n] + y_current[n+1]) - # stiff wave equation y_xx = y_current[n+1] - 2*y_current[n] + y_current[n-1] y_xxxx = y_current[n-2] - 4*y_current[n-1] + 6*y_current[n] - 4*y_current[n+1] + y_current[n+2] - term1 = (2 - 2*sigma*dt) * y_current[n] - term2 = (2*sigma*dt - 1) * y_last[n] - term3 = r2 * y_xx - s2 * y_xxxx - y_next[n] = term1 + term2 + term3 + y_next[n] = a1 * y_current[n] + a2 * y_last[n] + r2 * y_xx - s2 * y_xxxx y_next[0] = 0 y_next[1] = 0 @@ -101,8 +102,15 @@ for i in range(2, I): y_last = y_current.copy() y_current = y_next.copy() - if(i % 1000 == 0): - print(i/I * 100, "% complete") + if(i == 120000): + applyImpulse() + + if(i % 10000*seconds == 0): + print(f"{i/I * 100:4.4}% complete") + +end_time = time.perf_counter() +elapsed = end_time - start_time +print(f"Executed in {elapsed:.3f} seconds. {elapsed/seconds*100:.2f}% overshoot") plt.plot(np.arange(0, I, 1), y_sample) plt.grid() diff --git a/src/synth/Instrument.hpp b/src/synth/Instrument.hpp index 72249e3..f504e67 100644 --- a/src/synth/Instrument.hpp +++ b/src/synth/Instrument.hpp @@ -14,14 +14,14 @@ public: Instrument(ConfigService* config, LoggerService* logger); ~Instrument() = default; - void noteOn(float frequency, float velocity); - void noteOff(); + virtual void noteOn(float frequency, float velocity); + virtual void noteOff(); - bool isActive(); + virtual bool isActive(); - float process(bool& scopeTrigger); + virtual float process(bool& scopeTrigger); -private: +protected: float sampleRate_ = 44100.0f; bool active_ = false; diff --git a/src/synth/Instruments/PianoString.cpp b/src/synth/Instruments/PianoString.cpp new file mode 100644 index 0000000..28371ce --- /dev/null +++ b/src/synth/Instruments/PianoString.cpp @@ -0,0 +1,101 @@ + +#include "PianoString.hpp" + +PianoString::PianoString(ConfigService* config, LoggerService* logger) : Instrument(config, logger) { + + stringY_current_.resize(segmentCount_ + 1); + stringY_previous_.resize(segmentCount_ + 1); + stringY_next_.resize(segmentCount_ + 1); + stringX_.resize(segmentCount_ + 1); + +} + +void PianoString::noteOn(float frequency, float velocity) { + + logger_->log("Piano", LogFlag::Debug, "Note On"); + + // resize the state vectors so that they are stable at the specified frequency + stringY_current_.resize(segmentCount_ + 1); + stringY_previous_.resize(segmentCount_ + 1); + stringY_next_.resize(segmentCount_ + 1); + stringX_.resize(segmentCount_ + 1); + + // because stringGrid is resized, reevaluate + dx_ = stringLength_ / static_cast(segmentCount_); + for(size_t i = 0; i < segmentCount_ + 1; i++) { + stringX_[i] = i * dx_; + } + + // first iteration + for(size_t i = 2; i < segmentCount_ - 2; i++) { + stringY_current_[i] = 0.5f * r1_*r1_ * (stringY_previous_[i-1] - 2.0f*stringY_previous_[i] + stringY_previous_[i+1]); + } + // apply the velocity impulse + for(size_t i = 0; i < segmentCount_ + 1; i++) { + float v0 = impulseVelocity_ * std::exp(-1.0f * (stringX_[i] - strikePosition_)*(stringX_[i] - strikePosition_) / ((2.0f * impulseWidth_)*(2.0f * impulseWidth_))); + stringY_current_[i] = stringY_current_[i] + dt_ * v0; + } + + damping_ = 0.5f; + rms_ = 0.5f; + + // recalculate based on change in damping + a1_ = 2.0f - 2.0f * damping_ * dt_; + a2_ = 2.0f * damping_ * dt_ - 1.0f; +} + +void PianoString::noteOff() { + logger_->log("Piano", LogFlag::Debug, "Note Off"); + damping_ = 10.0f; + // recalculate based on change in damping + a1_ = 2.0f - 2.0f * damping_ * dt_; + a2_ = 2.0f * damping_ * dt_ - 1.0f; +} + +bool PianoString::isActive() { + return (std::abs(rms_) > 0.001f); +} + +float PianoString::process(bool& scopeTrigger) { + + /* + for n in range(2, N-2): + # stiff wave equation + y_xx = y_current[n+1] - 2*y_current[n] + y_current[n-1] + y_xxxx = y_current[n-2] - 4*y_current[n-1] + 6*y_current[n] - 4*y_current[n+1] + y_current[n+2] + + y_next[n] = a1 * y_current[n] + a2 * y_last[n] + r2 * y_xx - s2 * y_xxxx + + y_next[0] = 0 + y_next[1] = 0 + y_next[N-1] = 0 + y_next[N-2] = 0 + + # y_sample[i] = math.tanh(y_next[n_sample]) + y_sample[i] = y_next[n_sample] + + y_last = y_current.copy() + y_current = y_next.copy() + */ + + // simulate over string + for(size_t i = 2; i < segmentCount_ - 2; i++) { + float y_xx = stringY_current_[i-1] - 2.0f*stringY_current_[i] + stringY_current_[i+1]; + float y_xxxx = stringY_current_[i-2] - 4.0f*stringY_current_[i-1] + 6.0f*stringY_current_[i] - 4.0f*stringY_current_[i+1] + stringY_current_[i+2]; + stringY_next_[i] = a1_ * stringY_current_[i] + a2_ * stringY_previous_[i] + r2_ * y_xx - s2_ * y_xxxx; + } + stringY_next_[0] = 0.0f; + stringY_next_[1] = 0.0f; + stringY_next_[segmentCount_-1] = 0.0f; + stringY_next_[segmentCount_-2] = 0.0f; + + stringY_previous_ = stringY_current_; + stringY_current_ = stringY_next_; + + float sampleOut = stringY_next_[static_cast(samplePosition_*stringLength_*segmentCount_)]; + + rms_ = 0.99f * rms_ + 0.01f * sampleOut*sampleOut; + + return sampleOut; + +} \ No newline at end of file diff --git a/src/synth/Instruments/PianoString.hpp b/src/synth/Instruments/PianoString.hpp new file mode 100644 index 0000000..f693ee6 --- /dev/null +++ b/src/synth/Instruments/PianoString.hpp @@ -0,0 +1,68 @@ + +#pragma once + +#include + +#include "synth/Instrument.hpp" + +class PianoString : public Instrument { + +public: + + PianoString() = default; + PianoString(ConfigService* config, LoggerService* logger); + ~PianoString() = default; + + void noteOn(float frequency, float velocity) override; + void noteOff() override; + + bool isActive() override; + + float process(bool& scopeTrigger) override; + +private: + + // states + std::vector stringY_current_; + std::vector stringY_previous_; + std::vector stringY_next_; + std::vector stringX_; + + // constants + + // string parameters + size_t segmentCount_ = 30; + static constexpr float rho_ = 8000.0f; // density, steel, kg/m^3 + static constexpr float radius_ = 0.001f; // meters + static constexpr float stringTension_ = 1200.0f; // string tension, N + static constexpr float stiffness_ = 0.001f; // stiffness coefficient + float damping_ = 0.5f; // damping coefficient + static constexpr float stringLength_ = 1.0f; // length of string + static constexpr float strikePosition_ = 0.2f; // x of impulse location + static constexpr float impulseWidth_ = 0.02f; // x of impulse width + static constexpr float impulseVelocity_ = 10000.0f; // x/t of impulse magnitude + static constexpr float samplePosition_ = 0.1f; // percentage along L of sampling for audio + float crossSectionalArea_ = pi * std::pow(radius_, 2.0f); // string cross sectional area, assuming circular + float mu_ = crossSectionalArea_ * rho_; // linear mass density + float waveVelocity_ = std::sqrt(stringTension_ / mu_); // transverse wave velocity + + // eventually we'll have to dynamically tune our string according to the note that comes in + // an alternative is a fully built piano and then it calls voices under the instrument instead of how we do it currently + float f0_ = waveVelocity_ / (2.0f * stringLength_); // fundamental frequency of a non-stiff string + float f1_ = f0_ * std::sqrt(1.0f + stiffness_); // fundamental frequency of the stiff string + + float dx_ = stringLength_ / static_cast(segmentCount_); + float dt_ = 1.0f / sampleRate_; + + // derived constants + float r1_ = waveVelocity_ * dt_/dx_; + float r2_ = std::pow(waveVelocity_ * dt_/dx_, 2.0f); + float s1_ = stiffness_ * dt_/std::pow(dx_, 2.0f); + float s2_ = std::pow(stiffness_ * dt_/std::pow(dx_, 2.0f), 2.0f); + float a1_ = 2.0f - 2.0f * damping_ * dt_; + float a2_ = 2.0f * damping_ * dt_ - 1.0f; + + // keeping track of the string's activeness + float rms_ = 0.0f; + +}; diff --git a/src/synth/Voice.cpp b/src/synth/Voice.cpp index 6f96a5c..6b87ead 100644 --- a/src/synth/Voice.cpp +++ b/src/synth/Voice.cpp @@ -4,7 +4,8 @@ Voice::Voice(ConfigService* config, LoggerService* logger) : config_(config), logger_(logger) { - instrument_ = Instrument(config_, logger_); + // TODO: instrument factory + instrument_ = PianoString(config_, logger_); } diff --git a/src/synth/Voice.hpp b/src/synth/Voice.hpp index 9484580..f215050 100644 --- a/src/synth/Voice.hpp +++ b/src/synth/Voice.hpp @@ -4,6 +4,7 @@ #include #include "Instrument.hpp" +#include "Instruments/PianoString.hpp" // a voice is a tone generator that the synth uses for polyphony // the synth mixes multiple voices together into a polyphonic audio. calculations for samples are handled in the instrument @@ -37,6 +38,6 @@ private: ConfigService* config_; LoggerService* logger_; - Instrument instrument_; + PianoString instrument_; };