From ae2bac160139e93ccceed190aed1d404e055212e Mon Sep 17 00:00:00 2001 From: Gary Scavone Date: Tue, 16 Apr 2019 11:04:41 -0400 Subject: [PATCH] Added new Recorder class, fixed bug in StifKarp. --- include/Recorder.h | 167 +++++++++++++++ include/StifKarp.h | 3 +- projects/demo/Makefile.in | 4 +- projects/demo/tcl/Demo.tcl | 19 ++ projects/demo/tcl/Physical.tcl | 366 +++++++++++++++++---------------- projects/demo/utilities.cpp | 6 +- src/Flute.cpp | 4 +- src/Makefile.in | 2 +- src/Recorder.cpp | 345 +++++++++++++++++++++++++++++++ src/StifKarp.cpp | 4 +- 10 files changed, 736 insertions(+), 184 deletions(-) create mode 100644 include/Recorder.h create mode 100644 src/Recorder.cpp diff --git a/include/Recorder.h b/include/Recorder.h new file mode 100644 index 0000000..d9ca5a6 --- /dev/null +++ b/include/Recorder.h @@ -0,0 +1,167 @@ +#ifndef STK_RECORDER_H +#define STK_RECORDER_H + +#include "Instrmnt.h" +#include "Iir.h" +#include "DelayL.h" +#include "Noise.h" +#include "SineWave.h" +#include "ADSR.h" + +namespace stk { + +/***************************************************/ +/*! \class Recorder + \brief A recorder / flute physical model. + + This class implements a physical model of a recorder / + flute instrument, based on the paper "Sound production + in recorderlike instruments. II. A simulation model." + by M.P. Verge, A. Hirschberg and R. Causse, Journal of + the Acoustical Society of America, 1997. + + Control Change Numbers: + - Softness = 2 + - Noise Gain = 4 + - Noise Cutoff = 16 + - Vibrato Frequency = 11 + - Vibrato Gain = 1 + - Breath Pressure = 128 + + by Mathias Bredholt, McGill University. + Formatted for STK by Gary Scavone, 2019. +*/ +/***************************************************/ + +class Recorder : public Instrmnt +{ + public: + //! Class constructor. + Recorder( void ); + + //! Class destructor. + ~Recorder( void ); + + //! Reset and clear all internal state. + void clear( void ); + + //! Set instrument parameters for a particular frequency. + void setFrequency( StkFloat val ); + + //! Apply breath velocity to instrument with given amplitude and rate of increase. + void startBlowing( StkFloat amplitude, StkFloat rate ); + + //! Decrease breath velocity with given rate of decrease. + void stopBlowing( StkFloat rate ); + + //! Start a note with the given frequency and amplitude. + void noteOn( StkFloat frequency, StkFloat amplitude ); + + //! Stop a note with the given amplitude (speed of decay). + void noteOff( StkFloat amplitude ); + + //! Perform the control change specified by \e number and \e value (0.0 - 128.0). + void controlChange( int number, StkFloat value ); + + //! Compute and return one output sample. + StkFloat tick( unsigned int channel = 0 ); + + //! Fill a channel of the StkFrames object with computed outputs. + /*! + The \c channel argument must be less than the number of + channels in the StkFrames argument (the first channel is specified + by 0). However, range checking is only performed if _STK_DEBUG_ + is defined during compilation, in which case an out-of-range value + will trigger an StkError exception. + */ + StkFrames& tick( StkFrames& frames, unsigned int channel = 0 ); + + void setBlowPressure( StkFloat val ); + void setVibratoGain( StkFloat val ); + void setVibratoFrequency( StkFloat val ); + void setNoiseGain( StkFloat val ); + void setBreathCutoff( StkFloat val ); + void setSoftness( StkFloat val ); + + private: + DelayL pinDelay_; + DelayL poutDelay_; + DelayL jetDelay_; + Iir radiation_filter_; + Iir visco_in_filter_; + Iir visco_out_filter_; + Iir jetFilter_; + Noise turb_; + Iir turbFilter_; + SineWave vibrato_; + ADSR adsr_; + + //StkFloat M{ 0 }; + //StkFloat maxPressure_( 0 ); + double maxPressure_; + //StkFloat blow{ 0 }; + StkFloat vibratoGain_; + StkFloat noiseGain_; + StkFloat breathCutoff_; + StkFloat outputGain_; + StkFloat psi_; + + StkFloat poutL_; + StkFloat pout_; + StkFloat poutm1_; + StkFloat poutm2_; + StkFloat pin_; + StkFloat pinm1_; + StkFloat pinm2_; + + StkFloat b1; + StkFloat b3; + StkFloat b4; + + StkFloat Uj_; + StkFloat Ujm1_; + + StkFloat Qj_; + StkFloat Qjm1_; + StkFloat Qjm2_; + + StkFloat Q1_; + StkFloat Q1m1_; + StkFloat Q1m2_; + + StkFloat Qp_; + StkFloat Qpm1_; + + StkFloat pm_; +}; + +inline StkFrames& Recorder :: tick( StkFrames& frames, unsigned int channel ) +{ + unsigned int nChannels = lastFrame_.channels(); +#if defined(_STK_DEBUG_) + if ( channel > frames.channels() - nChannels ) { + oStream_ << "Recorder::tick(): channel and StkFrames arguments are incompatible!"; + handleError( StkError::FUNCTION_ARGUMENT ); + } +#endif + + StkFloat *samples = &frames[channel]; + unsigned int j, hop = frames.channels() - nChannels; + if ( nChannels == 1 ) { + for ( unsigned int i=0; i b_coeffs( &b_rad[0], &b_rad[0]+3 ); + std::vector a_coeffs( &a_rad[0], &a_rad[0]+3 ); + radiation_filter_.setCoefficients(b_coeffs, a_coeffs); + + // Visco-thermal loss filter + StkFloat b_visco[4] = { 0.83820223947141, -0.16888603248373, -0.64759781930259, 0.07424498608506 }; + StkFloat a_visco[4] = { 1.0, -0.33623476246554, -0.71257915055968, 0.14508304017256 }; + b_coeffs.clear(); + b_coeffs.assign( &b_visco[0], &b_visco[0]+4 ); + a_coeffs.clear(); + a_coeffs.assign( &a_visco[0], &a_visco[0]+4 ); + visco_in_filter_.setCoefficients(b_coeffs, a_coeffs); + visco_out_filter_.setCoefficients(b_coeffs, a_coeffs); + + setBreathCutoff( 500 ); + setFrequency( 880 ); +} + +Recorder :: ~Recorder( void ) +{ +} + +void Recorder :: clear( void ) +{ + pinDelay_.clear(); + poutDelay_.clear(); + jetDelay_.clear(); + radiation_filter_.clear(); + visco_in_filter_.clear(); + visco_out_filter_.clear(); + turbFilter_.clear(); +} + +void Recorder :: setFrequency( StkFloat val ) +{ +#if defined(_STK_DEBUG_) + if ( val <= 0.0 ) { + oStream_ << "Recorder::setFrequency: argument is less than or equal to zero!"; + handleError( StkError::WARNING ); return; + } +#endif + + StkFloat M = Stk::sampleRate() / val - 4 - 3; + pinDelay_.setDelay( M ); + poutDelay_.setDelay( M ); +} + +void Recorder :: setBlowPressure( StkFloat val ) +{ + maxPressure_ = val; +} + +void Recorder :: setVibratoGain( StkFloat val ) +{ + vibratoGain_ = val; +} + +void Recorder :: setVibratoFrequency( StkFloat val ) +{ + vibrato_.setFrequency( val ); +} + +void Recorder :: setNoiseGain( StkFloat val ) +{ + noiseGain_ = val; +} + +void Recorder :: setBreathCutoff( StkFloat val ) +{ + // The gain of this filter is quite high + breathCutoff_ = val; + StkFloat Q = 0.99; + StkFloat r = 2.0 * sin(PI * val / sampleRate()); + StkFloat q = 1.0 - r * Q; + StkFloat as[3] = { 1.0, r * r - q - 1, q }; + std::vector b_turb(1, r*r); + std::vector a_turb( &as[0], &as[0]+3 ); + turbFilter_.setCoefficients(b_turb, a_turb); +} + +void Recorder :: setSoftness( StkFloat val ) +{ + psi_ = val; +} + +void Recorder :: startBlowing( StkFloat amplitude, StkFloat rate ) +{ + if ( amplitude <= 0.0 || rate <= 0.0 ) { + oStream_ << "Recorder::startBlowing: one or more arguments is less than or equal to zero!"; + handleError( StkError::WARNING ); return; + } + + adsr_.setAttackRate( rate ); + //maxPressure_ = amplitude / (StkFloat) 0.8; + maxPressure_ = 35 * amplitude; + adsr_.keyOn(); +} + +void Recorder :: stopBlowing( StkFloat rate ) +{ + if ( rate <= 0.0 ) { + oStream_ << "Recorder::stopBlowing: argument is less than or equal to zero!"; + handleError( StkError::WARNING ); return; + } + + adsr_.setReleaseRate( rate ); + adsr_.keyOff(); +} + +void Recorder :: noteOn( StkFloat frequency, StkFloat amplitude ) +{ + this->setFrequency( frequency ); + this->startBlowing( 1.1 + (amplitude * 0.20), amplitude * 0.02 ); + outputGain_ = amplitude / 40.0; +} + +void Recorder :: noteOff( StkFloat amplitude ) +{ + this->stopBlowing( amplitude * 0.02 ); +} + +void Recorder :: controlChange( int number, StkFloat value ) +{ +#if defined(_STK_DEBUG_) + if ( Stk::inRange( value, 0.0, 128.0 ) == false ) { + oStream_ << "Recorder::controlChange: value (" << value << ") is out of range!"; + handleError( StkError::WARNING ); return; + } +#endif + + StkFloat normalizedValue = value * ONE_OVER_128; + if (number == 2) // 2 + psi_ = 2.0 * normalizedValue; + else if (number == 16) + setBreathCutoff( normalizedValue * 2000 ); + else if (number == __SK_NoiseLevel_) // 4 + noiseGain_ = normalizedValue; + else if (number == __SK_ModFrequency_) // 11 + vibrato_.setFrequency( normalizedValue * 12.0); + else if (number == __SK_ModWheel_) // 1 + vibratoGain_ = ( normalizedValue * 0.4 ); + else if (number == __SK_AfterTouch_Cont_) // 128 + maxPressure_ = 35.0 * normalizedValue; +#if defined(_STK_DEBUG_) + else { + oStream_ << "Recorder::controlChange: undefined control number (" << number << ")!"; + handleError( StkError::WARNING ); + } +#endif +} + +StkFloat Recorder::tick( unsigned int ) +{ + // Read in from delay lines + pinm2_ = pinm1_; + pinm1_ = pin_; + pin_ = pinDelay_.lastOut(); + + poutm2_ = poutm1_; + poutm1_ = pout_; + poutL_ = poutDelay_.lastOut(); + + // Filter wave components for visco-thermal losses + pin_ = visco_in_filter_.tick(pin_); + poutL_ = visco_out_filter_.tick(poutL_); + + // Get input blow pressure + StkFloat pf = maxPressure_ * adsr_.tick() * (vibrato_.tick() * vibratoGain_ + (1 - vibratoGain_)); + + StkFloat T = 1.0 / sampleRate(); + + // Jet velocity at flue exit + Ujm1_ = Uj_; + Uj_ = Ujm1_ + T / (rho * lc) * (pf - pm_ - 0.5 * rho * Ujm1_ * Ujm1_); + + // Jet flow at flue exit + Qjm2_ = Qjm1_; + Qjm1_ = Qj_; + Qj_ = h * H * Uj_; + + // Jet drive + StkFloat Uj_steady = fmax(sqrt(2 * pf / rho), 0.1); + StkFloat fc_jet = 0.36 / W * Uj_steady; + StkFloat g_jet = 0.002004 * exp(-0.06046 * Uj_steady); + StkFloat r_jet = 0.95 - Uj_steady * 0.015; + StkFloat b0_jet = g_jet * (1 - r_jet * r_jet) / 2; + + // Calculate coefficients for resonant filter + StkFloat b_jet[3] = { b0_jet, 0, -b0_jet }; + StkFloat a_jet[3] = { 1, -2 * r_jet * cos(2 * PI * fc_jet * T), r_jet * r_jet }; + std::vector b_jetcoeffs( &b_jet[0], &b_jet[0]+3 ); + std::vector a_jetcoeffs( &a_jet[0], &a_jet[0]+3 ); + jetFilter_.setCoefficients( b_jetcoeffs, a_jetcoeffs ); + StkFloat eta = jetFilter_.tick(jetDelay_.lastOut()); + + // Calculate flow source Q1 + Q1m1_ = Q1_; + Q1_ = b * H * Uj_ * (1 + tanh(eta / (psi_ * b))); + + // Calculate pressure pulse modeling the jet drive + StkFloat pjd = -rho * dd / Sm * (Q1_ - Q1m1_) / T; + + // Vortex shedding + int Qp_sign = 0; + if (Qp_ < 0) Qp_sign = -1; + else if (Qp_ > 0) Qp_sign = 1; + + StkFloat pa = -0.5 * rho * (Qp_ / (0.6 * Sm)) * (Qp_ / (0.6 * Sm)) * Qp_sign; + + // Turbulence + StkFloat pt = turbFilter_.tick(noiseGain_ * turb_.tick() * 0.5 * rho * Uj_ * Uj_); + + // Pressure pulse delta p + StkFloat dp = pjd + pa + pt; + + // Calculate outgoing pressure pout + pout_ = ((b3 - b1 * b2 - 1) * pin_ + + (2 * b1 * b2 - b3) * (pinm1_ - poutm1_) + + b1 * b2 * (poutm2_ - pinm2_) - + b1 * (Qj_ - 2 * Qjm1_ + Qjm2_) + + b4 * (Qj_ - Qjm1_) + dp) / (1 - b1 * b2 + b3); + + // Flow in the pipe + Qpm1_ = Qp_; + Qp_ = Sp / (rho * c0) * (pout_ - pin_); + + // Mouth pressure + pm_ = pout_ + pin_ - dp + rho * din / Sm * (Qp_ - Qpm1_)/T; + + // Calculate transverse acoustic velocity + StkFloat Q1d = Q1_ - 0.5 * b * H * Uj_; + StkFloat Vac = 2.0 / PI * Qp_ / Sm - 0.38 * Q1d / Sm; + jetDelay_.tick(Vac); + + // Calculate new jet delay line length + //jet_.setDelay(fmin(W / (0.6 * Uj_steady) * sampleRate(), 200.0)); + jetDelay_.setDelay(fmin(W / (0.6 * Uj_steady * T), 200.0)); + + // Radiation loss filtering + StkFloat pin_L = radiation_filter_.tick(poutL_); + + // Write to delay lines + poutDelay_.tick(pout_); + pinDelay_.tick(pin_L); + + lastFrame_[0] = outputGain_ * (pout_ + pin_); + return lastFrame_[0]; + + //return (pout_0 + pin_0) * 0.01; +} + +} // stk namespace diff --git a/src/StifKarp.cpp b/src/StifKarp.cpp index 2275f6a..b9358f5 100644 --- a/src/StifKarp.cpp +++ b/src/StifKarp.cpp @@ -89,7 +89,7 @@ void StifKarp :: setStretch( StkFloat stretch ) StkFloat freq = lastFrequency_ * 2.0; StkFloat dFreq = ( (0.5 * Stk::sampleRate()) - freq ) * 0.25; StkFloat temp = 0.5 + (stretch * 0.5); - if ( temp > 0.9999 ) temp = 0.9999; + if ( temp > 0.99999 ) temp = 0.99999; for ( int i=0; i<4; i++ ) { coefficient = temp * temp; biquad_[i].setA2( coefficient ); @@ -131,7 +131,7 @@ void StifKarp :: pluck( StkFloat amplitude ) } pluckAmplitude_ = amplitude; - for ( unsigned long i=0; itick() * pluckAmplitude_) );