oscillatorsEOC.lib

// =============================================================================
// ========== oscillatorsEOC.lib ===============================================
// =============================================================================
//
// This library contains band-limited oscillators with arbitrary harmonic 
// content for classic analogue waveforms such as sawtooth, square, 
// and triangle waves, as well as band-limited pulse-trains with arbitrary
// duty-cycles. The library includes a quadrature oscillator based on,
// arguably, the best recursive design available in the literature.
//
// The environment prefix is "o2".
//
// List of functions:
//
//    blit_bi,
//    blit_bi_duty,
//    blit_uni,
//    osc_quad,
//    ph,
//    pulse_train,
//    saw,
//    square,
//    tri.
//
// Copyright (c) 2019-2020, Dario Sanfilippo <sanfilippo.dario at gmail dot com>
// All rights reserved.

declare name "Oscillators Library";
declare author "Dario Sanfilippo";
declare copyright "Copyright (c) 2019-2020, Dario Sanfilippo <sanfilippo.dario
      at gmail dot com>";
declare version "2.0.0";
declare license "GPLv2.0";

au = library("auxiliaryEOC.lib");
ba = library("basics.lib");
fi = library("filters.lib");
f2 = library("filtersEOC.lib");
ma = library("maths.lib");
m2 = library("mathsEOC.lib");
os = library("oscillators.lib");
o2 = library("oscillatorsEOC.lib");
ro = library("routes.lib");
si = library("signals.lib");

// o2.blit_bi(H[n], F[n]); -----------------------------------------------------
//
// Bipolar band-limited impulse train (BLIT) based on period sinc function. 
//
// The maximum harmonic number is given by:
//
//    rint(SR/frequency/4).
//
// The bipolar BLIT has no DC component and the number of harmonics include 
// the fundamental frequency and its odd multiples. Variations in the harmonic 
// content take place at the beginning of each cycle to avoid clicks. Hence, 
// changes will take place after a time that is the period of the BLIT or less.
//
// Unlike the technique described in [Stilson and Smith 1996] where the 
// bipolar BLIT is implemented by summation of a unipolar BLIT with its 
// delayed and inverted copy, the technique showed here uses an even ratio 
// between the frequencies of the sine functions used to generate the sinc, 
// which results in a correct harmonic content (odd harmonics) for any given 
// BLIT frequencies.
//
// The amplitude of the function, regardless of the harmonics, 
// is normalised to unit-amplitude peaks.
//
// 2 inputs:
//    H[n], (rounded to rint inside the asinc function) number of harmonics;
//    F[n], BLIT frequency in Hz.
//
// 1 outputs:
//    y[n], bipolar band-limited impulse train.
//
blit_bi(h, f) = m2.asinc_bi(h1, phase)
      with {
           lim = rint(m2.div(ma.SR, f) / 4);
           h1 = ba.sAndH(trigger, min(lim, h));
           trigger = (ma.signum(f) * (phase - phase') < 0);
           // Add au.dirac to "trigger" if you require the initial value
           // of H[n] to be triggered initially at n = 0.
           phase = os.phasor(1, f);
      };
// -----------------------------------------------------------------------------

// o2.blit_bi_duty(H[n], D[n], F[n]); ------------------------------------------
//
// Bipolar band-limited impulse train with arbitrary duty cycle
// following the paper by [Stilson and Smith 1996]:
//
// https://ccrma.stanford.edu/~stilti/papers/blit.pdf.
//
// The lowest frequency at which an entire duty cycle can be explored is
// 1 Hz. If lower frequencies are required, the first argument of
// m2.delta should be changed.
//
// 3 inputs:
//    H[n], (rounded to rint) number of harmonics;
//    D[n], duty cycle in the range [0; 1];
//    F[n], frequency of the BLIT in Hz.
//
// 1 outputs:
//    y[n], bipolar BLIT with arbitrary duty cycle.
//
blit_bi_duty(h, d, f) = m2.delta(1, d1, o2.blit_uni(h, f))
      with {
           d1 = d * m2.div(1, f);
      };
// -----------------------------------------------------------------------------

// o2.blit_uni(H[n], F[n]); ----------------------------------------------------
//
// The unipolar BLIT has a DC component and the number of harmonics include the 
// fundamental frequency and its multiples (both even and odd). 
//
// The technique described here is based on the paper by
// [Stilson and Smith 1996]:
//
// https://ccrma.stanford.edu/~stilti/papers/blit.pdf.
//
// The amplitude of the function, regardless of the harmonics, 
// is normalised to unit-amplitude peaks.
//
// 2 inputs:
//    H[n], (rounded to rint inside the asinc function) number of harmonics;
//    F[n], BLIT frequency in Hz.
//
// 1 outputs:
//    y[n], unipolar band-limited impulse train.
//
blit_uni(h, f) = m2.asinc_uni(h1, phase)
      with {
           lim = floor(m2.div(ma.SR, f) / 2);
           h1 = ba.sAndH(trigger, min(lim, h));
           trigger = (ma.signum(f) * (phase - phase') < 0);
           // Add au.dirac to "trigger" if you require the initial value
           // of H[n] to be triggered initially at n = 0.
           phase = os.phasor(1, f);
      };
// -----------------------------------------------------------------------------

// o2.osc_quad(F[n]); ----------------------------------------------------------
//
// Recursive quadrature oscillator by Martin Vicanek. This design is
// arguably the best recursive quadrature oscillator available in the
// literature. The system shows long-term stability as well as accuracy at
// low frequencies. 
//
// Ref: https://vicanek.de/articles/QuadOsc.pdf.
//
// 1 inputs:
//    F[n], oscillator frequency in Hz.
//
// 2 outputs:
//    y1[n], cosine (real part);
//    y2[n], sine (imaginary part).
//
osc_quad(f) =  tick 
               ~ ( _ , 
                   _)
      with {
           k1 = tan(ma.PI * f / ma.SR);
           k2 = 2 * k1 / (1 + k1 * k1);
           tick(u, v) =    omega - k1 * (v + k2 * omega) , 
                           v + k2 * omega
           with {
               omega = (u + au.dirac) - k1 * v;
           };
      };
// -----------------------------------------------------------------------------

// o2.ph(F[n], R[n]); ----------------------------------------------------------
//
// Phasor with reset input.
//
// Note: the arguments of the function should be inverted for
// consistency, but that requires checking for backward compatibility.
//
// 2 inputs:
//    F[n], frequency in Hz;
//    R[n], reset phasor to zero if R[n] != 0.
//
// 1 outputs:
//    y[n], phasor output.
//
ph(freq, reset) =  (+ (freq / ma.SR * r) : ma.decimal)
                   ~ * (r)
      with {
           r = reset == 0;
      };
// -----------------------------------------------------------------------------

// o2.pulse_train(H[n], D[n], F[n]); -------------------------------------------
//
// BLIT-based variable width pulse train. Implemented following
// [Stilson and Smith 1996]:
//
// https://ccrma.stanford.edu/~stilti/papers/blit.pdf.
//
// 3 inputs:
//    H[n], number of harmonics, both even and odd, cast to the closest
//             INT and phase-locked to the beginning of each cycle. The harmonic
//             content is affected by the duty cycle;
//    D[n], duty cycle in the range [0; 1];
//    F[n], pulse train frequency in Hz.
//
// 1 outputs:
//    y[n], band-limited pulse train.
//
pulse_train(h, d, f) = 
      blit_bi_duty(h, d, f) : m2.div(f2.leaky(.1 / m2.twopi), scale) + d
      with {
           lim = rint(m2.div(ma.SR, f) / 2);
           scale = m2.div(lim, min(lim, h)) : si.smooth(ba.tau2pole(.1));
      };
// -----------------------------------------------------------------------------

// o2.saw(H[n], F[n]); ---------------------------------------------------------
//
// BLIT-based band-limited sawtooth oscillator.
//
// This function is based on [Stilson and Smith 1996], although the
// scaling factor for a unit-amplitude normalisation has been determined by 
// the author.
//
// Ref: https://ccrma.stanford.edu/~stilti/papers/blit.pdf.
//
// 2 inputs:
//    H[n], number of harmonics (both even and odd), cast to the closest INT;
//    F[n], frequency of the oscillator in Hz.
//
// 1 outputs:
//    y[n], band-limited sawtooth oscillator.
//
saw(h, f) = blit_uni(h, f) <: _ - f2.lp1p(f / 100) : 
      m2.div(f2.leaky(.1 / m2.twopi), scale) : fi.highpass(1, 20)
      with {
           lim = floor(m2.div(ma.SR, f) / 2);
           scale = m2.div(lim, (2 * min(lim, h))) : si.smooth(ba.tau2pole(.1));
      };
// -----------------------------------------------------------------------------

// o2.square(H[n], F[n]); ------------------------------------------------------
//
// BLIT-based band-limited square oscillator. This technique implements
// the integration of a bipolar BLIT. The bipolar BLIT is based on a sinc
// function with even ratios between the sine functions used in sinc, which
// result in a more precise harmonic content throughout the entire
// frequency range. The scaling factor normalises the output to unit
// amplitude for all frequencies and harmonics.
//
// 2 inputs:
//    H[n], number of harmonics (odd), cast to the closest
//             INT and phase-locked to the beginning of each cycle.
//    F[n], frequency of the oscillator in Hz.
//
// 1 outputs:
//    y[n], band-limited square oscillator.
//
square(h, f) = blit_bi(h, f) : m2.div(f2.leaky(.1 / m2.twopi), scale) : 
      fi.highpass(1, 20)
      with {
           lim = rint(m2.div(ma.SR, f) / 4);
           scale = m2.div(lim, (2 * min(lim, h))) : si.smooth(ba.tau2pole(.1));
      };
// -----------------------------------------------------------------------------

// o2.triangle(H[n], F[n]); ----------------------------------------------------
//
// BLIT-based band-limited triangle oscillator. This technique implements
// the integration of a BLIT-based square wave. The bipolar BLIT in the square 
// wave is based on a sinc function with even ratios between the sine functions 
// used in sinc, which result in a more precise harmonic content throughout 
// the entire frequency range. The scaling factor normalises the output to unit
// amplitude for all frequencies and harmonics
//
// 2 inputs:
//    H[n], number of harmonics (odd), cast to the closest
//             INT and phase-locked to the beginning of each cycle.
//    F[n], frequency of the oscillator in Hz.
//
// 1 outputs:
//    y[n], band-limited triangle oscillator.
//
triangle(h, f) = square(h, f) : m2.div(f2.leaky(.1 / m2.twopi), scale) : 
      fi.highpass(1, 20)
      with {
           scale = rint(m2.div(ma.SR, f) / 4) : si.smooth(ba.tau2pole(.1));
      };
// -----------------------------------------------------------------------------