Added new Recorder class, fixed bug in StifKarp.

src/Flute.cpp

@@ -47,8 +47,8 @@ Flute :: Flute( StkFloat lowestFrequency )
   adsr_.setAllTimes( 0.005, 0.01, 0.8, 0.010 );
   endReflection_ = 0.5;
   jetReflection_ = 0.5;
-  noiseGain_     = 0.15;    // Breath pressure random component.
-  vibratoGain_   = 0.05;    // Breath periodic vibrato component.
+  noiseGain_     = 0.15;    // Breath pressure random component
+  vibratoGain_   = 0.05;    // Breath periodic vibrato component
   jetRatio_      = 0.32;
 
 	maxPressure_ = 0.0;

src/Makefile.in

@@ -26,7 +26,7 @@ OBJECTS	=	Stk.o Generator.o Noise.o Blit.o BlitSaw.o BlitSquare.o Granulate.o \
 					Function.o ReedTable.o JetTable.o BowTable.o Cubic.o \
 					Voicer.o Vector3D.o Sphere.o Twang.o Guitar.o \
 					\
-					Instrmnt.o Clarinet.o BlowHole.o Saxofony.o Flute.o Brass.o BlowBotl.o \
+					Instrmnt.o Clarinet.o BlowHole.o Saxofony.o Flute.o Recorder.o Brass.o BlowBotl.o \
 					Bowed.o Plucked.o StifKarp.o Sitar.o Mandolin.o Mesh2D.o \
 					FM.o Rhodey.o Wurley.o TubeBell.o HevyMetl.o PercFlut.o BeeThree.o FMVoices.o \
 					Sampler.o Moog.o Simple.o Drummer.o Shakers.o \

src/Recorder.cpp

@@ -0,0 +1,345 @@
+/***************************************************/
+/*! \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.
+*/
+/***************************************************/
+
+#include "Recorder.h"
+#include "SKINImsg.h"
+
+namespace stk {
+
+// Air constants
+const stk::StkFloat rho = 1.2041; // density of air
+const stk::StkFloat c0 = 343.21;  // speed of sound in air
+
+// Flute constants
+const stk::StkFloat lc = 0.02;     // length of flue canal
+const stk::StkFloat h = 0.001;     // height of flue exit
+const stk::StkFloat H = 0.02;      // pipe diameter
+const stk::StkFloat W = 4 * h;     // width of mouth
+const stk::StkFloat Sp = H * H;    // cross-section of pipe
+const stk::StkFloat Sm = W * H;    // cross-section of mouth
+const stk::StkFloat din = 0.0030;  // end correction
+const stk::StkFloat dout = 0.0063; // end correction
+const stk::StkFloat dm = din + dout;  // end correction of mouth
+const stk::StkFloat dd = 0.0035;   // acoustic distance between Q1 and Q2
+const stk::StkFloat rp = sqrt(Sp / stk::PI);
+const stk::StkFloat b = 0.4 * h;   // jet width
+
+// Calculation coefficients
+const stk::StkFloat b2 = Sp / (rho * c0);
+
+Recorder :: Recorder()
+{
+  vibratoGain_ = 0.0;
+  noiseGain_ = 0.2;
+  breathCutoff_ = 0.0;
+  outputGain_ = 0.0;
+  psi_ = 1.0;
+  poutL_ = 0;
+  pout_ = 0;
+  poutm1_ = 0;
+  poutm2_ = 0;
+  pin_ = 0;
+  pinm1_ = 0;
+  pinm2_ = 0;
+  Uj_ = 0;
+  Ujm1_ = 0;
+  Qj_ = 0;
+  Qjm1_ = 0;
+  Qjm2_ = 0;
+  Q1_ = 0;
+  Q1m1_ = 0;
+  Q1m2_ = 0;
+  Qp_ = 0;
+  Qpm1_ = 0;
+  pm_ = 0;
+  
+  pinDelay_.tick( 0 );
+  poutDelay_.tick( 0 );
+  jetDelay_.tick( 0 );
+
+  jetDelay_.setDelay(200);
+  vibrato_.setFrequency(4);
+
+  // Calculation coefficients ... would need to be recalculated if sample rate changes
+  StkFloat T = 1.0 / Stk::sampleRate();
+  b1 = rho / (4.0 * PI * c0 * T * T);
+  b3 = dm * Sp / (T * Sm * c0);
+  b4 = rho * dout / (Sm * T);
+
+  // Radiation loss filter
+  StkFloat A = rp * rp / (4 * c0 * c0 * T * T);
+  StkFloat B = 0.82 * rp / (c0*T);
+  StkFloat b_rad[3] = { 1 + A - B, B - 2 * A, A };
+  StkFloat a_rad[3] = { A - B - 1, B - 2 * A, A };
+  std::vector<StkFloat> b_coeffs( &b_rad[0], &b_rad[0]+3 );
+  std::vector<StkFloat> 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<StkFloat> b_turb(1, r*r);
+  std::vector<StkFloat> 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<StkFloat> b_jetcoeffs( &b_jet[0], &b_jet[0]+3 );
+  std::vector<StkFloat> 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

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; i<length_; i++ ) {
+  for ( unsigned long i=0; i<lastLength_; i++ ) {
     // Fill delay with noise additively with current contents.
     delayLine_.tick( (delayLine_.lastOut() * 0.6) + 0.4 * noise_.tick() * pluckAmplitude_ );
     //delayLine_.tick( combDelay_.tick((delayLine_.lastOut() * 0.6) + 0.4 * noise->tick() * pluckAmplitude_) );