Add BLEP-corrected hard sync

Source/Synthesis/oscillator.cpp

@@ -20,8 +20,11 @@ float Oscillator::Process()
         case WAVE_POLYBLEP_TRI:
             t   = phase_;
             out = phase_ < 0.5f ? 1.0f : -1.0f;
-            out += Polyblep(phase_inc_, t);
-            out -= Polyblep(phase_inc_, fastmod1f(t + 0.5f));
+            if(!skip_builtin_polyblep_)
+            {
+                out += Polyblep(phase_inc_, t);
+                out -= Polyblep(phase_inc_, fastmod1f(t + 0.5f));
+            }
             // Leaky Integrator:
             // y[n] = A + x[n] + (1 - A) * y[n-1]
             out       = phase_inc_ * out + (1.0f - phase_inc_) * last_out_;
@@ -31,23 +34,39 @@ float Oscillator::Process()
         case WAVE_POLYBLEP_SAW:
             t   = phase_;
             out = (2.0f * t) - 1.0f;
-            out -= Polyblep(phase_inc_, t);
+            if(!skip_builtin_polyblep_)
+                out -= Polyblep(phase_inc_, t);
             out *= -1.0f;
             break;
         case WAVE_POLYBLEP_SQUARE:
             t   = phase_;
             out = phase_ < pw_ ? 1.0f : -1.0f;
-            out += Polyblep(phase_inc_, t);
-            out -= Polyblep(phase_inc_, fastmod1f(t + (1.0f - pw_)));
+            if(!skip_builtin_polyblep_)
+            {
+                out += Polyblep(phase_inc_, t);
+                out -= Polyblep(phase_inc_, fastmod1f(t + (1.0f - pw_)));
+            }
             out *= 0.707f; // ?
             break;
         default: out = 0.0f; break;
     }
+
+    // Apply any pending sync-BLEP residual (set by SyncReset(frac) on the
+    // previous sample). Consumed in one shot.
+    out += sync_corr_;
+    sync_corr_             = 0.0f;
+    skip_builtin_polyblep_ = false;
+
     phase_ += phase_inc_;
     if(phase_ > 1.0f)
     {
         phase_ -= 1.0f;
         eoc_ = true;
+        // Fractional position within this sample at which the wrap occurred.
+        // Derivation: before advance, phase was p_before = phase_ + 1 - inc.
+        // Wrap time f satisfies p_before + f*inc = 1, so
+        //   f = (1 - p_before) / inc = 1 - phase_/inc.
+        eoc_frac_ = (phase_inc_ > 0.0f) ? (1.0f - phase_ / phase_inc_) : 0.0f;
     }
     else
     {
@@ -58,6 +77,63 @@ float Oscillator::Process()
     return out * amp_;
 }
 
+float Oscillator::NaiveWaveformValue(float phase) const
+{
+    float t;
+    switch(waveform_)
+    {
+        case WAVE_SIN: return sinf(phase * TWOPI_F);
+        case WAVE_TRI:
+            t = -1.0f + (2.0f * phase);
+            return 2.0f * (fabsf(t) - 0.5f);
+        case WAVE_SAW: return -1.0f * ((phase * 2.0f) - 1.0f);
+        case WAVE_RAMP: return (phase * 2.0f) - 1.0f;
+        case WAVE_SQUARE: return phase < pw_ ? 1.0f : -1.0f;
+        // For the polyblep variants we return the *pre-correction* naive value
+        // in the same sign/scale convention as Process()'s output (minus the
+        // polyblep residual and, for TRI, minus the leaky integrator).
+        case WAVE_POLYBLEP_TRI: return phase < 0.5f ? 1.0f : -1.0f;
+        case WAVE_POLYBLEP_SAW: return -1.0f * ((2.0f * phase) - 1.0f);
+        case WAVE_POLYBLEP_SQUARE:
+            return (phase < pw_ ? 1.0f : -1.0f) * 0.707f;
+        default: return 0.0f;
+    }
+}
+
+void Oscillator::SyncReset(float frac)
+{
+    // Clamp frac into [0, 1) defensively.
+    if(frac < 0.0f) frac = 0.0f;
+    if(frac > 1.0f) frac = 1.0f;
+
+    // Slave's pre-sync continuous phase at the instant sync fired.
+    float pre_phase = phase_ + frac * phase_inc_;
+    if(pre_phase >= 1.0f) pre_phase -= 1.0f;
+    float pre_val = NaiveWaveformValue(pre_phase);
+
+    // Slave's phase at the next sample boundary: it reset to 0 at the sync
+    // instant and has been running for (1 - frac) * phase_inc_ since.
+    float post_phase = (1.0f - frac) * phase_inc_;
+    float post_val   = NaiveWaveformValue(post_phase);
+
+    // Step magnitude at the discontinuity.
+    float step = post_val - pre_val;
+
+    // Polyblep residual for a unit step at fractional offset `frac` into the
+    // period [k-1, k], evaluated at sample k (the sample right after the
+    // step). Standard 2-sample polyblep shape gives correction = -frac^2/2
+    // per unit of step magnitude at sample k. The pre-step sample k-1 has
+    // already been emitted, so its correction is lost (causality limit of
+    // an online implementation).
+    sync_corr_             = step * (-0.5f * frac * frac);
+    skip_builtin_polyblep_ = true;
+
+    phase_    = post_phase;
+    last_out_ = 0.0f;
+    eoc_      = true;
+    eoc_frac_ = 0.0f;
+}
+
 float Oscillator::CalcPhaseInc(float f)
 {
     return f * sr_recip_;

Source/Synthesis/oscillator.h

@@ -49,16 +49,19 @@ class Oscillator
     */
     void Init(float sample_rate)
     {
-        sr_        = sample_rate;
-        sr_recip_  = 1.0f / sample_rate;
-        freq_      = 100.0f;
-        amp_       = 0.5f;
-        pw_        = 0.5f;
-        phase_     = 0.0f;
-        phase_inc_ = CalcPhaseInc(freq_);
-        waveform_  = WAVE_SIN;
-        eoc_       = true;
-        eor_       = true;
+        sr_                     = sample_rate;
+        sr_recip_               = 1.0f / sample_rate;
+        freq_                   = 100.0f;
+        amp_                    = 0.5f;
+        pw_                     = 0.5f;
+        phase_                  = 0.0f;
+        phase_inc_              = CalcPhaseInc(freq_);
+        waveform_               = WAVE_SIN;
+        eoc_                    = true;
+        eor_                    = true;
+        eoc_frac_               = 0.0f;
+        sync_corr_              = 0.0f;
+        skip_builtin_polyblep_  = false;
     }
 
 
@@ -92,6 +95,14 @@ class Oscillator
     */
     inline bool IsEOC() { return eoc_; }
 
+    /** When IsEOC() is true, returns the fractional position within the
+        just-processed sample at which the phase wrap occurred.
+        0.0 = at the very start of the sample period, 1.0 = at the end.
+        Used by a slave oscillator's SyncReset(frac) for sub-sample-accurate
+        (BLEP-corrected) hard sync.
+    */
+    inline float GetEocFraction() const { return eoc_frac_; }
+
     /** Returns true if cycle rising.
     */
     inline bool IsRising() { return phase_ < 0.5f; }
@@ -112,13 +123,31 @@ class Oscillator
     */
     void Reset(float _phase = 0.0f) { phase_ = _phase; }
 
+    /** BLEP-corrected hard sync. `frac` is the sub-sample position (0..1)
+        at which the master wrapped, typically obtained from
+        `master.GetEocFraction()`. Produces a band-limited reset: the slave's
+        phase is advanced as if it had been reset at the correct sub-sample
+        instant, and a polyblep residual is stored to be applied to the next
+        Process() output. Suppresses the built-in polyblep on that sample
+        since its assumed step magnitude would be incorrect for a sync event.
+
+        Intended usage:
+            float m = master.Process();
+            if(master.IsEOC()) slave.SyncReset(master.GetEocFraction());
+            float s = slave.Process();
+    */
+    void SyncReset(float frac);
+
   private:
     float   CalcPhaseInc(float f);
+    float   NaiveWaveformValue(float phase) const;
     uint8_t waveform_;
     float   amp_, freq_, pw_;
     float   sr_, sr_recip_, phase_, phase_inc_;
     float   last_out_, last_freq_;
+    float   eoc_frac_, sync_corr_;
     bool    eor_, eoc_;
+    bool    skip_builtin_polyblep_;
 };
 } // namespace daisysp
 #endif