//############################## oscillators.lib ######################################
// This library contains a collection of sound generators. Its official prefix is `os`.
//
// The oscillators library is organized into 9 sections:
//
// * [Wave-Table-Based Oscillators](#wave-table-based-oscillators)
// * [Low Frequency Oscillators](#low-frequency-oscillators)
// * [Low Frequency Sawtooths](#low-frequency-sawtooths)
// * [Alias-Suppressed Sawtooth](#alias-suppressed-sawtooth)
// * [Alias-Suppressed Pulse, Square, and Impulse Trains](#alias-suppressed-pulse-square-and-impulse-trains)
// * [Filter-Based Oscillators](#filter-based-oscillators)
// * [Waveguide-Resonator-Based Oscillators](#waveguide-resonator-based-oscillators)
// * [Casio CZ Oscillators](#casio-cz-oscillators)
// * [PolyBLEP-Based Oscillators](#polyblep-based-oscillators)
//
// #### References
// * <https://github.com/grame-cncm/faustlibraries/blob/master/oscillators.lib>
//########################################################################################

/************************************************************************
************************************************************************
FAUST library file, GRAME section

Except where noted otherwise, Copyright (C) 2003-2017 by GRAME,
Centre National de Creation Musicale.
----------------------------------------------------------------------
GRAME LICENSE

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation; either version 2.1 of the
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA.

EXCEPTION TO THE LGPL LICENSE : As a special exception, you may create a
larger FAUST program which directly or indirectly imports this library
file and still distribute the compiled code generated by the FAUST
compiler, or a modified version of this compiled code, under your own
copyright and license. This EXCEPTION TO THE LGPL LICENSE explicitly
grants you the right to freely choose the license for the resulting
compiled code. In particular the resulting compiled code has no obligation
to be LGPL or GPL. For example you are free to choose a commercial or
closed source license or any other license if you decide so.
************************************************************************
************************************************************************/

ma = library("maths.lib");
ba = library("basics.lib");
fi = library("filters.lib");
si = library("signals.lib");

declare name "Faust Oscillator Library";
declare version "1.6.0";

// This library contains platform specific constants
pl = library("platform.lib");

//======================Oscillators based on mathematical functions===============
//
// Note that there is a numerical problem with several phasor functions built using the internal
// `phasor_imp`. The reason is that the incremental step is smaller than `ma.EPSILON`, which happens with very small frequencies, 
// so it will have no effect when summed to 1, but it will be enough to make the fractional function wrap 
// around when summed to 0. An example of this problem can be observed when running the following code:
//
// `process = os.phasor(1.0, -.001);`
//
// The output of this program is the sequence 1, 0, 1, 0, 1... This happens because the negative incremental 
// step is greater than `-ma.EPSILON`, which will have no effect when summed to 1, but it will be significant 
// enough to make the fractional function  wrap around when summed to 0.
//
// The incremental step can be clipped to guarantee that the phasor will 
// always run correctly for its full cycle, otherwise, for increments smaller than `ma.EPSILON`, 
// phasor would initially run but it'd eventually get stuck once the output gets big enough.
// 
// All functions using `phasor_imp` are affected by this problem, but a safer
// version is implemented, and can be used alternatively by setting `SAFE=1` in the environment using 
// [explicit sustitution](https://faustdoc.grame.fr/manual/syntax/#explicit-substitution) syntax.
// 
// For example: `process = os[SAFE=1;].phasor(1.0, -.001);` will use the safer implementation of `phasor_imp`.
//=================================================================================

//=========================Wave-Table-Based Oscillators===================================
// Oscillators using tables. The table size is set by the 
// [pl.tablesize](https://github.com/grame-cncm/faustlibraries/blob/master/platform.lib) constant.
//========================================================================================

// Global parameter to use the safer version of `phasor_imp`, but which
// could be used in other functions as well.

SAFE = 0; // 0: use the faster version, 1: use the safer version

//-----------------------`(os.)sinwaveform`------------------------
// Sine waveform ready to use with a `rdtable`.
//
// #### Usage
//
// ```
// sinwaveform(tablesize) : _
// ```
//
// Where:
//
// * `tablesize`: the table size
//------------------------------------------------------------
sinwaveform(tablesize) =
    sin(float(ba.period(tablesize)) * (2.0 * ma.PI) / float(tablesize));


//-----------------------`(os.)coswaveform`------------------------
// Cosine waveform ready to use with a `rdtable`.
//
// #### Usage
//
// ```
// coswaveform(tablesize) : _
// ```
//
// Where:
//
// * `tablesize`: the table size
//------------------------------------------------------------
coswaveform(tablesize) =
    cos(float(ba.period(tablesize)) * (2.0 * ma.PI) / float(tablesize));

// Possibly faster version using integer arithmetic
phasor_env(freq, N) = environment {

    //------- GLOBAL PARAMS
    nbits = 31;
    tablesize = 1<<N;
    accuracy = int(nbits - N);
    mask = (1<<nbits)-1;
    inc(step) = int(tablesize * step * (1<<accuracy));

    //------- LAMBDA DSP CASE 
    lambda(inc_op) = (inc_op : &(mask)) ~ _ : >>(accuracy) : /(tablesize);

    //------- MINIMAL CASE
    hsp(0,0) = lambda(+(inc(freq/ma.SR)'));
    
    //------- GENERAL CASE
    hsp(reset,phase) = lambda(select2(hard_reset,+(inc(freq/ma.SR)),inc(phase)))
    with {
        hard_reset = (1-1')|reset;
    };
};
declare phasor_env author "Pierre Mascarade Relano, Maxime Sirbu, Stéphane Letz";

// Generic phasor with `reset` and `phase` parameters to be specialised in concrete use-cases.
phasor_imp(freq, reset, phase) = (select2(hard_reset, +(incr(SAFE)), phase) : ma.decimal) ~ _
with {
    incr_aux = freq/ma.SR;

    // Faster but less accurate version
    incr(0) = incr_aux;

    // To make sure that the incremental step is greater or equal to EPSILON or 
    // less than or equal to -EPSILON to avoid numerical problems. 
    // A frequency of 0Hz can still be used to freeze the phasor.
    incr(1)= (freq != 0) * ba.if(freq < 0, min(-1.0 * ma.EPSILON, incr_aux), max(ma.EPSILON, incr_aux));

    // To correctly start at `phase` at the first sample
    hard_reset = (1-1')|reset; 
};

// Possibly faster version using integer arithmetic
// phasor_imp(freq, reset, phase) = phasor_env(freq, 16).hsp(reset, phase);

// Version to be used with tables
phasor_table(tablesize, freq, reset, phase) = phasor_imp(freq, reset, phase) : *(float(tablesize));


//-----------------------`(os.)phasor`------------------------
// A simple phasor to be used with a `rdtable`.
// `phasor` is a standard Faust function.
//
// #### Usage
//
// ```
// phasor(tablesize,freq) : _
// ```
//
// Where:
//
// * `tablesize`: the table size
// * `freq`: the frequency in Hz
//
// Note that `tablesize` is just a multiplier for the output of a unit-amp phasor
// so `phasor(1.0, freq)` can be used to generate a phasor output in the range [0, 1[.
//------------------------------------------------------------
phasor(tablesize, freq) = phasor_table(tablesize, freq, 0, 0);


//-----------------------`(os.)hs_phasor`------------------------
// Hardsyncing phasor to be used with a `rdtable`.
//
// #### Usage
//
// ```
// hs_phasor(tablesize,freq,reset) :  _
// ```
//
// Where:
//
// * `tablesize`: the table size
// * `freq`: the frequency in Hz
// * `reset`: a reset signal, reset phase to 0 when equal to 1
//---------------------------------------------------------
declare hs_phasor author "Mike Olsen, revised by Stéphane Letz";

hs_phasor(tablesize, freq, reset) = phasor_table(tablesize, freq, reset, 0);


//-----------------------`(os.)hsp_phasor`------------------------
// Hardsyncing phasor with selectable phase to be used with a `rdtable`.
//
// #### Usage
//
// ```
// hsp_phasor(tablesize,freq,reset,phase)
// ```
//
// Where:
//
// * `tablesize`: the table size
// * `freq`: the frequency in Hz
// * `reset`: reset the oscillator to phase when equal to 1
// * `phase`: phase between 0 and 1
//---------------------------------------------------------
declare hsp_phasor author "Christophe Lebreton, revised by Stéphane Letz";

hsp_phasor(tablesize, freq, reset, phase) = phasor_table(tablesize, freq, reset, phase);


//-----------------------`(os.)oscsin`------------------------
// Sine wave oscillator.
// `oscsin` is a standard Faust function.
//
// #### Usage
//
// ```
// oscsin(freq) : _
// ```
//
// Where:
//
// * `freq`: the frequency in Hz
//------------------------------------------------------------
oscsin(freq) = rdtable(tablesize, sinwaveform(tablesize), int(phasor(tablesize,freq)))
with {
    tablesize = pl.tablesize;
};


//-----------------------`(os.)hs_oscsin`------------------------
// Sin lookup table with hardsyncing phase.
//
// #### Usage
//
// ```
// hs_oscsin(freq,reset) : _
// ```
//
// Where:
//
// * `freq`: the frequency in Hz
// * `reset`: reset the oscillator to 0 when equal to 1
//---------------------------------------------------------
declare hs_oscsin author "Mike Olsen";

hs_oscsin(freq,reset) = rdtable(tablesize, sinwaveform(tablesize), int(hs_phasor(tablesize,freq,reset)))
with {
    tablesize = pl.tablesize;
};


//-----------------------`(os.)osccos`------------------------
// Cosine wave oscillator.
//
// #### Usage
//
// ```
// osccos(freq) : _
// ```
//
// Where:
//
// * `freq`: the frequency in Hz
//------------------------------------------------------------
osccos(freq) = rdtable(tablesize, coswaveform(tablesize), int(phasor(tablesize,freq)))
with {
    tablesize = pl.tablesize;
};


//-----------------------`(os.)hs_osccos`------------------------
// Cos lookup table with hardsyncing phase.
//
// #### Usage
//
// ```
// hs_osccos(freq,reset) : _
// ```
//
// Where:
//
// * `freq`: the frequency in Hz
// * `reset`: reset the oscillator to 0 when equal to 1
//---------------------------------------------------------
declare hs_osccos author "Stéphane Letz";

hs_osccos(freq,reset) = rdtable(tablesize, coswaveform(tablesize), int(hs_phasor(tablesize,freq,reset)))
with {
    tablesize = pl.tablesize;
};


//-----------------------`(os.)oscp`------------------------
// A sine wave generator with controllable phase.
//
// #### Usage
//
// ```
// oscp(freq,phase) : _
// ```
//
// Where:
//
// * `freq`: the frequency in Hz
// * `phase`: the phase in radian
//------------------------------------------------------------
oscp(freq,phase) = oscsin(freq) * cos(phase) + osccos(freq) * sin(phase);


//-----------------------`(os.)osci`------------------------
// Interpolated phase sine wave oscillator.
//
// #### Usage
//
// ```
// osci(freq) : _
// ```
//
// Where:
//
// * `freq`: the frequency in Hz
//------------------------------------------------------------
osci(freq) = s1 + d * (s2 - s1)
with {
    tablesize = pl.tablesize;
    i = int(phasor(tablesize,freq));
    d = ma.decimal(phasor(tablesize,freq));
    s1 = rdtable(tablesize+1,sinwaveform(tablesize),i);
    s2 = rdtable(tablesize+1,sinwaveform(tablesize),i+1);
};


//-----------------------`(os.)osc`------------------------
// Default sine wave oscillator (same as [oscsin](#oscsin)).
// `osc` is a standard Faust function.
//
// #### Usage
//
// ```
// osc(freq) : _
// ```
//
// Where:
//
// * `freq`: the frequency in Hz
//------------------------------------------------------------
osc = oscsin;


//-----------------------`(os.)m_oscsin`------------------------
// Sine wave oscillator based on the `sin` mathematical function.
//
// #### Usage
//
// ```
// m_oscsin(freq) : _
// ```
//
// Where:
//
// * `freq`: the frequency in Hz
//------------------------------------------------------------
m_oscsin(freq) = lf_sawpos(freq) : *(2*ma.PI) : sin;


//-----------------------`(os.)m_osccos`------------------------
// Sine wave oscillator based on the `cos` mathematical function.
//
// #### Usage
//
// ```
// m_osccos(freq) : _
// ```
//
// Where:
//
// * `freq`: the frequency in Hz
//------------------------------------------------------------
m_osccos(freq) = lf_sawpos(freq) : *(2*ma.PI) : cos;


// end GRAME section
//########################################################################################
/************************************************************************
FAUST library file, jos section

Except where noted otherwise, The Faust functions below in this
section are Copyright (C) 2003-2022 by Julius O. Smith III <jos@ccrma.stanford.edu>
([jos](http://ccrma.stanford.edu/~jos/)), and released under the
(MIT-style) [STK-4.3](#stk-4.3-license) license.

The MarkDown comments in this section are Copyright 2016-2017 by Romain
Michon and Julius O. Smith III, and are released under the
[CCA4I](https://creativecommons.org/licenses/by/4.0/) license (TODO: if/when Romain agrees)

************************************************************************/

//===============================Low Frequency Oscillators===============================
// Low Frequency Oscillators (LFOs) have prefix `lf_`
// (no aliasing suppression, since it is inaudible at LF).
// Use `sawN` and its derivatives for audio oscillators with suppressed aliasing.
//==================================================================

//--------`(os.)lf_imptrain`----------
// Unit-amplitude low-frequency impulse train.
// `lf_imptrain` is a standard Faust function.

// #### Usage
//
// ```
// lf_imptrain(freq) : _
// ```

// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
lf_imptrain(freq) = lf_sawpos(freq)<:-(mem)<0; // definition below


//--------`(os.)lf_pulsetrainpos`----------
// Unit-amplitude nonnegative LF pulse train, duty cycle between 0 and 1.
//
//
// #### Usage
//
// ```
// lf_pulsetrainpos(freq, duty) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
// * `duty`: duty cycle between 0 and 1
//------------------------------------------------------------
lf_pulsetrainpos(freq,duty) = float(lf_sawpos(freq) <= duty);

//pulsetrainpos = lf_pulsetrainpos; // for backward compatibility


//--------`(os.)lf_pulsetrain`----------
// Unit-amplitude zero-mean LF pulse train, duty cycle between 0 and 1.
//
// #### Usage
//
// ```
// lf_pulsetrain(freq,duty) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
// * `duty`: duty cycle between 0 and 1
//------------------------------------------------------------
lf_pulsetrain(freq,duty) = 2.0*lf_pulsetrainpos(freq,duty) - 1.0;


//--------`(os.)lf_squarewavepos`----------
// Positive LF square wave in [0,1]
//
// #### Usage
//
// ```
// lf_squarewavepos(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
lf_squarewavepos(freq) = lf_pulsetrainpos(freq,0.5);
// squarewavepos = lf_squarewavepos; // for backward compatibility


//--------`(os.)lf_squarewave`----------
// Zero-mean unit-amplitude LF square wave.
// `lf_squarewave` is a standard Faust function.
//
// #### Usage
//
// ```
// lf_squarewave(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
lf_squarewave(freq) = 2.0*lf_squarewavepos(freq) - 1.0;
// squarewave = lf_squarewave; // for backward compatibility


//--------`(os.)lf_trianglepos`----------
// Positive unit-amplitude LF positive triangle wave.
//
// #### Usage
//
// ```
// lf_trianglepos(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
lf_trianglepos(freq) = 1.0-abs(saw1(freq)); // saw1 defined below


//----------`(os.)lf_triangle`----------
// Zero-mean unit-amplitude LF triangle wave.
// `lf_triangle` is a standard Faust function.
//
// #### Usage
//
// ```
// lf_triangle(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
declare lf_triangle author "Bart Brouns";
declare lf_triangle licence "STK-4.3";

lf_triangle(freq) = 2.0*lf_trianglepos(freq) - 1.0;


//================== Low Frequency Sawtooths ====================
// Sawtooth waveform oscillators for virtual analog synthesis et al.
// The 'simple' versions (`lf_rawsaw`, `lf_sawpos` and `saw1`), are mere samplings of
// the ideal continuous-time ("analog") waveforms.  While simple, the
// aliasing due to sampling is quite audible.  The differentiated
// polynomial waveform family (`saw2`, `sawN`, and derived functions)
// do some extra processing to suppress aliasing (not audible for
// very low fundamental frequencies).  According to Lehtonen et al.
// (JASA 2012), the aliasing of `saw2` should be inaudible at fundamental
// frequencies below 2 kHz or so, for a 44.1 kHz sampling rate and 60 dB SPL
// presentation level;  fundamentals 415 and below required no aliasing
// suppression (i.e., `saw1` is ok).
//=====================================================================

//-----------------`(os.)lf_rawsaw`--------------------
// Simple sawtooth waveform oscillator between 0 and period in samples.
//
// #### Usage
//
// ```
// lf_rawsaw(periodsamps) : _
// ```
//
// Where:
//
// * `periodsamps`: number of periods per samples
//---------------------------------------------------------
lf_rawsaw(periodsamps) = (_,periodsamps : fmod) ~ +(1.0);


//-----------------`(os.)lf_sawpos`--------------------
// Simple sawtooth waveform oscillator between 0 and 1.
//
// #### Usage
//
// ```
// lf_sawpos(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//
//---------------------------------------------------------
declare lf_sawpos author "Bart Brouns, revised by Stéphane Letz";
declare lf_sawpos licence "STK-4.3";

lf_sawpos(freq) = phasor_imp(freq, 0, 0);


//-----------------`(os.)lf_sawpos_phase`--------------------
// Simple sawtooth waveform oscillator between 0 and 1
// with phase control.
//
// #### Usage
//
// ```
// lf_sawpos_phase(freq, phase) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
// * `phase`: phase between 0 and 1
//---------------------------------------------------------
declare lf_sawpos_phase author "Bart Brouns, revised by Stéphane Letz";
declare lf_sawpos_phase licence "STK-4.3";

lf_sawpos_phase(freq,phase) = phasor_imp(freq, 0, phase);


//-----------------`(os.)lf_sawpos_reset`--------------------
// Simple sawtooth waveform oscillator between 0 and 1
// with reset.
//
// #### Usage
//
// ```
// lf_sawpos_reset(freq,reset) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
// * `reset`: reset the oscillator to 0 when equal to 1
//
//---------------------------------------------------------
declare lf_sawpos_reset author "Bart Brouns, revised by Stéphane Letz";
declare lf_sawpos_reset licence "STK-4.3";

lf_sawpos_reset(freq,reset) = phasor_imp(freq, reset, 0);


//-----------------`(os.)lf_sawpos_phase_reset`--------------------
// Simple sawtooth waveform oscillator between 0 and 1
// with phase control and reset.
//
// #### Usage
//
// ```
// lf_sawpos_phase_reset(freq,phase,reset) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
// * `phase`: phase between 0 and 1
// * `reset`: reset the oscillator to phase when equal to 1
//
//---------------------------------------------------------
declare lf_sawpos_phase_reset author "Bart Brouns, revised by Stéphane Letz";
declare lf_sawpos_phase_reset licence "STK-4.3";

lf_sawpos_phase_reset(freq,phase,reset) = phasor_imp(freq, reset, phase);


//-----------------`(os.)lf_saw`--------------------
// Simple sawtooth waveform oscillator between -1 and 1.
// `lf_saw` is a standard Faust function.
//
// #### Usage
//
// ```
// lf_saw(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//---------------------------------------------------------
declare saw1 author "Bart Brouns";
declare saw1 licence "STK-4.3";

saw1(freq) = 2.0 * lf_sawpos(freq) - 1.0;
lf_saw(freq) = saw1(freq);

//================== Alias-Suppressed Sawtooth ====================
//-----------------`(os.)sawN`--------------------
// Alias-Suppressed Sawtooth Audio-Frequency Oscillator using Nth-order polynomial transitions
// to reduce aliasing.
//
// `sawN(N,freq)`, `sawNp(N,freq,phase)`, `saw2dpw(freq)`, `saw2(freq)`, `saw3(freq)`,
// `saw4(freq)`, `sawtooth(freq)`, `saw2f2(freq)`, `saw2f4(freq)`
//
// #### Usage
//
// ```
// sawN(N,freq) : _        // Nth-order aliasing-suppressed sawtooth using DPW method (see below)
// sawNp(N,freq,phase) : _ // sawN with phase offset feature
// saw2dpw(freq) : _       // saw2 using DPW
// saw2ptr(freq) : _       // saw2 using the faster, stateless PTR method
// saw2(freq) : _          // DPW method, but subject to change if a better method emerges
// saw3(freq) : _          // sawN(3)
// saw4(freq) : _          // sawN(4)
// sawtooth(freq) : _      // saw2
// saw2f2(freq) : _        // saw2dpw with 2nd-order droop-correction filtering
// saw2f4(freq) : _        // saw2dpw with 4th-order droop-correction filtering
// ```
//
// Where:
//
// * `N`: polynomial order, a constant numerical expression between 1 and 4
// * `freq`: frequency in Hz
// * `phase`: phase between 0 and 1
//
// #### Method
// Differentiated Polynomial Wave (DPW).
//
// ##### Reference
// "Alias-Suppressed Oscillators based on Differentiated Polynomial Waveforms",
// Vesa Valimaki, Juhan Nam, Julius Smith, and Jonathan Abel,
// IEEE Tr. Audio, Speech, and Language Processing (IEEE-ASLP),
// Vol. 18, no. 5, pp 786-798, May 2010.
// 10.1109/TASL.2009.2026507.
//
// #### Notes
// The polynomial order `N` is limited to 4 because noise has been
// observed at very low `freq` values.  (LFO sawtooths should of course
// be generated using `lf_sawpos` instead.)
//-----------------------------------------------------------------
declare sawN author "Julius O. Smith III";
declare sawN license "STK-4.3";
// --- sawN for N = 1 to 4 ---
// Orders 5 and 6 have noise at low fundamentals: MAX_SAW_ORDER = 6; MAX_SAW_ORDER_NEXTPOW2 = 8;
MAX_SAW_ORDER = 4;
MAX_SAW_ORDER_NEXTPOW2 = 8; // par cannot handle the case of 0 elements
sawN(N,freq) = saw1l : poly(Nc) : D(Nc-1) : gate(Nc-1)
with {
  Nc = max(1,min(N,MAX_SAW_ORDER));
  clippedFreq = max(20.0,abs(freq)); // use lf_sawpos(freq) for LFOs (freq < 20 Hz)
  saw1l = 2*lf_sawpos(clippedFreq) - 1; // zero-mean, amplitude +/- 1
  poly(1,x) = x;
  poly(2,x) = x*x;
  poly(3,x) = x*x*x - x;
  poly(4,x) = x*x*(x*x - 2.0);
  poly(5,x) = x*(7.0/3 + x*x*(-10.0/3.0 + x*x));
  poly(6,x) = x*x*(7.0 + x*x*(-5.0 + x*x));
  p0n = float(ma.SR)/clippedFreq; // period in samples
  diff1(x) = (x - x')/(2.0/p0n);
  diff(N) = seq(n,N,diff1); // N diff1s in series
  factorial(0) = 1;
  factorial(i) = i * factorial(i-1);
  D(0) = _;
  D(i) = diff(i)/factorial(i+1);
  gate(N) = *(1@(N)); // delayed step for blanking startup glitch
};

//------------------`(os.)sawNp`--------------------------------
// Same as `(os.)sawN` but with a controllable waveform phase.
//
// #### Usage
//
// ```
// sawNp(N,freq,phase) : _
// ```
//
// where
//
// * `N`: waveform interpolation polynomial order 1 to 4 (constant integer expression)
// * `freq`: frequency in Hz
// * `phase`: waveform phase as a fraction of one period (rounded to nearest sample)
//
// #### Implementation Notes
//
// The phase offset is implemented by delaying `sawN(N,freq)` by
// `round(phase*ma.SR/freq)` samples, for up to 8191 samples.
// The minimum sawtooth frequency that can be delayed a whole period
// is therefore `ma.SR/8191`, which is well below audibility for normal
// audio sampling rates.
//
//-----------------------------------------------------------------
declare sawNp author "Julius O. Smith III";
declare sawNp license "STK-4.3";
// --- sawNp for N = 1 to 4 ---
// Phase offset = delay (max 8191 samples is more than one period of audio):
sawNp(N,freq,phase) = sawN(N,freq) : @(max(0,min(8191,int(0.5+phase*ma.SR/freq))));

//------------------`(os.)saw2, (os.)saw3, (os.)saw4`--------------
// Alias-Suppressed Sawtooth Audio-Frequency Oscillators of order 2, 3, 4.
//
// #### Usage
//
// ```
// saw2(freq) : _
// saw3(freq) : _
// saw4(freq) : _
// ```
//
// where
//
// * `freq`: frequency in Hz
//
// ##### References
// See `sawN` above.
//
// #### Implementation Notes
//
// Presently, only `saw2` uses the PTR method, while `saw3` and `saw4` use DPW.
// This is because PTR has been implemented and tested for the 2nd-order case only.
//
//------------------------------------------------------------------
saw2 = saw2ptr; // "faustlibraries choice"
saw3 = sawN(3); // only choice available right now
saw4 = sawN(4); // only choice available right now

//---------------------------`(os.)saw2ptr`---------------------------
// Alias-Suppressed Sawtooth Audio-Frequency Oscillator
// using Polynomial Transition Regions (PTR) for order 2.
//
// #### Usage
//
// ```
// saw2ptr(freq) : _
// ```
//
// where
//
// * `freq`: frequency in Hz
//
// ##### Implementation
//
// Polynomial Transition Regions (PTR) method for aliasing suppression.
//
// ##### References
//
// * Kleimola, J.; Valimaki, V., "Reducing Aliasing from Synthetic Audio
//     Signals Using Polynomial Transition Regions," in Signal Processing
//     Letters, IEEE , vol.19, no.2, pp.67-70, Feb. 2012
// * <https://aaltodoc.aalto.fi/bitstream/handle/123456789/7747/publication6.pdf?sequence=9>
// * <http://research.spa.aalto.fi/publications/papers/spl-ptr/>
//
// ##### Notes
//
// Method PTR may be preferred because it requires less
// computation and is stateless which means that the frequency `freq`
// can be modulated arbitrarily fast over time without filtering
// artifacts.  For this reason, `saw2` is presently defined as `saw2ptr`.
//
//--------------------------------------------------------
declare saw2ptr author "Julius O. Smith III";
declare saw2ptr license "STK-4.3";
// specialized reimplementation:
saw2ptr(freq) = y with { // newer PTR version (stateless - freq can vary at any speed)
  p0 = float(ma.SR)/float(max(ma.EPSILON,abs(freq))); // period in samples
  t0 = 1.0/p0; // phase increment
  p = ((_<:(-(1)<:_,_),_) <: selector1,selector2) ~(+(t0)):!,_;
  selector1 = select2(<(0)); // for feedback
  selector2 = select2(<(0), (_<:_,(*(1-p0):+(1)):+), _); // for output
  y = 2*p-1;
};

//----------------------`(os.)saw2dpw`---------------------
// Alias-Suppressed Sawtooth Audio-Frequency Oscillator
// using the Differentiated Polynomial Waveform (DWP) method.
//
// #### Usage
//
// ```
// saw2dpw(freq) : _
// ```
//
// where
//
// * `freq`: frequency in Hz
//
// This is the original Faust `saw2` function using the DPW method.
// Since `saw2` is now defined as `saw2ptr`, the DPW version
// is now available as `saw2dwp`.
//--------------------------------------------------------
declare saw2dpw author "Julius O. Smith III";
declare saw2dpw license "STK-4.3";
saw2dpw(freq) = saw1(freq) <: * <: -(mem) : *(0.25'*ma.SR/freq);

//------------------`(os.)sawtooth`--------------------------------
// Alias-suppressed aliasing-suppressed sawtooth oscillator, presently defined as `saw2`.
// `sawtooth` is a standard Faust function.
//
// #### Usage
//
// ```
// sawtooth(freq) : _
// ```
//
// with
//
// * `freq`: frequency in Hz
//--------------------------------------------------------
sawtooth = saw2; // default choice for sawtooth signal - see also sawN

//------------------`(os.)saw2f2, (os.)saw2f4`--------------------------------
// Alias-Suppressed Sawtooth Audio-Frequency Oscillator with Order 2 or 4 Droop Correction Filtering.
//
// #### Usage
//
// ```
// saw2f2(freq) : _
// saw2f4(freq) : _
// ```
//
// with
//
// * `freq`: frequency in Hz
//
// In return for aliasing suppression, there is some attenuation near half the sampling rate.
// This can be considered as beneficial, or it can be compensated with a high-frequency boost.
// The boost filter is second-order for `saw2f2` and fourth-order for `saw2f4`, and both are designed
// for the DWP case and therefore use `saw2dpw`.
// See Figure 4(b) in the DPW reference for a plot of the slight droop in the DPW case.
//--------------------------------------------------------
declare saw2f2 author "Julius O. Smith III";
declare saw2f2 license "STK-4.3";
// --- Correction-filtered versions of saw2: saw2f2, saw2f4 -----
saw2f2 = saw2dpw : cf2 with {
  cf2 = fi.tf2(1.155704605878911, 0.745184288225518,0.040305967265900,
        0.823765146386639, 0.117420665547108);
};
declare saw2f4 author "Julius O. Smith III";
declare saw2f4 license "STK-4.3";
saw2f4 = saw2dpw : cf4 with {
  cf4 = fi.iir((1.155727435125014, 2.285861038554662,
        1.430915027294021, 0.290713280893317, 0.008306401748854),
        (2.156834679164532, 1.559532244409321, 0.423036498118354,
        0.032080681130972));
};

//=========Alias-Suppressed Pulse, Square, and Impulse Trains============
// Alias-Suppressed Pulse, Square and Impulse Trains.
//
// `pulsetrainN`, `pulsetrain`, `squareN`, `square`, `imptrainN`, `imptrain`,
// `triangleN`, `triangle`
//
// All are zero-mean and meant to oscillate in the audio frequency range.
// Use simpler sample-rounded `lf_*` versions above for LFOs.
//
// #### Usage
//
// ```
// pulsetrainN(N,freq,duty) : _
// pulsetrain(freq, duty) : _ // = pulsetrainN(2)
//
// squareN(N,freq) : _
// square : _ // = squareN(2)
//
// imptrainN(N,freq) : _
// imptrain : _ // = imptrainN(2)
//
// triangleN(N,freq) : _
// triangle : _ // = triangleN(2)
// ```
//
// Where:
//
// * `N`: polynomial order, a constant numerical expression
// * `freq`: frequency in Hz
//====================================================================


//------------------`(os.)impulse`--------------------------------
// One-time impulse generated when the Faust process is started.
// `impulse` is a standard Faust function.
//
// #### Usage
//
// ```
// impulse : _
// ```
//--------------------------------------------------------
impulse = 1-1';


//------------------`(os.)pulsetrainN`--------------------------------
// Alias-suppressed pulse train oscillator.
//
// #### Usage
//
// ```
// pulsetrainN(N,freq,duty) : _
// ```
//
// Where:
//
// * `N`: order, as a constant numerical expression
// * `freq`: frequency in Hz
// * `duty`: duty cycle between 0 and 1

//--------------------------------------------------------
pulsetrainN(N,freq,duty) = diffdel(sawN(N,freqC),del) with {
 // non-interpolated-delay version: diffdel(x,del) = x - x@int(del+0.5);
 // linearly interpolated delay version (sounds good to me):
    diffdel(x,del) = x-x@int(del)*(1-ma.frac(del))-x@(int(del)+1)*ma.frac(del);
 // Third-order Lagrange interpolated-delay version (see filters.lib):
 // diffdel(x,del) = x - fdelay3(DELPWR2,max(1,min(DELPWR2-2,ddel)));
 DELPWR2 = 2048; // Needs to be a power of 2 when fdelay*() used above.
 delmax = DELPWR2-1; // arbitrary upper limit on diff delay (duty=0.5)
 SRmax = 96000.0; // assumed upper limit on sampling rate
 fmin = SRmax / float(2.0*delmax); // 23.4 Hz (audio freqs only)
 freqC = max(freq,fmin); // clip frequency at lower limit
 period = (float(ma.SR) / freqC); // actual period
 ddel = duty * period; // desired delay
 del = max(0,min(delmax,ddel));
};


//------------------`(os.)pulsetrain`--------------------------------
// Alias-suppressed pulse train oscillator. Based on `pulsetrainN(2)`.
// `pulsetrain` is a standard Faust function.
//
// #### Usage
//
// ```
// pulsetrain(freq,duty) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
// * `duty`: duty cycle between 0 and 1
//--------------------------------------------------------
pulsetrain = pulsetrainN(2);


//------------------`(os.)squareN`--------------------------------
// Alias-suppressed square wave oscillator.
//
// #### Usage
//
// ```
// squareN(N,freq) : _
// ```
//
// Where:
//
// * `N`: order, as a constant numerical expression
// * `freq`: frequency in Hz
//--------------------------------------------------------
squareN(N,freq) = pulsetrainN(N,freq,0.5);


//------------------`(os.)square`--------------------------------
// Alias-suppressed square wave oscillator. Based on `squareN(2)`.
// `square` is a standard Faust function.
//
// #### Usage
//
// ```
// square(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//--------------------------------------------------------
square = squareN(2);


//------------------`(os.)imptrainN`--------------------------------
// Alias-suppressed impulse train generator.
//
// #### Usage
//
// ```
// imptrainN(N,freq) : _
// ```
//
// Where:
//
// * `N`: order, as a constant numerical expression
// * `freq`: frequency in Hz
//--------------------------------------------------------
imptrainN(N,freq) = impulse + 0.5*ma.diffn(sawN(N,freq));


//------------------`(os.)imptrain`--------------------------------
// Alias-suppressed impulse train generator. Based on `imptrainN(2)`.
// `imptrain` is a standard Faust function.
//
// #### Usage
//
// ```
// imptrain(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//--------------------------------------------------------
imptrain = imptrainN(2); // default based on saw2


//------------------`(os.)triangleN`--------------------------------
// Alias-suppressed triangle wave oscillator.
//
// #### Usage
//
// ```
// triangleN(N,freq) : _
// ```
//
// Where:
//
// * `N`: order, as a constant numerical expression
// * `freq`: frequency in Hz
//--------------------------------------------------------
triangleN(N,freq) = squareN(N,freq) : fi.pole(p) : *(gain) with {
  gain = 4.0*freq/ma.SR; // for aproximate unit peak amplitude
  p = 0.999;
};


//------------------`(os.)triangle`--------------------------------
// Alias-suppressed triangle wave oscillator. Based on `triangleN(2)`.
// `triangle` is a standard Faust function.
//
// #### Usage
//
// ```
// triangle(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//--------------------------------------------------------
triangle = triangleN(2); // default based on saw2


//===============================Filter-Based Oscillators=================================
// Filter-Based Oscillators.
//
// #### Usage
//
// ```
// osc[b|rq|rs|rc|s](freq), where freq = frequency in Hz.
// ```
//
// #### References
//
// * <http://lac.linuxaudio.org/2012/download/lac12-slides-jos.pdf>
// * <https://ccrma.stanford.edu/~jos/pdf/lac12-paper-jos.pdf>
//========================================================================================

//--------------------------`(os.)oscb`--------------------------------
// Sinusoidal oscillator based on the biquad.
//
// #### Usage
//
// ```
// oscb(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
oscb(f) = impulse : fi.tf2(1,0,0,a1,1)
with {
  a1 = -2*cos(2*ma.PI*f/ma.SR);
};


//--------------------------`(os.)oscrq`---------------------------
// Sinusoidal (sine and cosine) oscillator based on 2D vector rotation,
//  = undamped "coupled-form" resonator
//  = lossless 2nd-order normalized ladder filter.
//
// #### Usage
//
// ```
// oscrq(freq) : _,_
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html>
//------------------------------------------------------------
oscrq(f) = impulse : fi.nlf2(f,1); // sine and cosine outputs

//--------------------------`(os.)oscrs`---------------------------
// Sinusoidal (sine) oscillator based on 2D vector rotation,
//  = undamped "coupled-form" resonator
//  = lossless 2nd-order normalized ladder filter.
//
// #### Usage
//
// ```
// oscrs(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html>
//------------------------------------------------------------
oscrs(f) = impulse : fi.nlf2(f,1) : _,!; // sine

//--------------------------`(os.)oscrc`---------------------------
// Sinusoidal (cosine) oscillator based on 2D vector rotation,
//  = undamped "coupled-form" resonator
//  = lossless 2nd-order normalized ladder filter.
//
// #### Usage
//
// ```
// oscrc(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html>
//------------------------------------------------------------
oscrc(f) = impulse : fi.nlf2(f,1) : !,_; // cosine

oscrp(f,p) = oscrq(f) : *(cos(p)), *(sin(p)) : + ; // p=0 for sine, p=PI/2 for cosine, etc.

oscr = oscrs; // default = sine (starts without a pop)

//--------------------------`(os.)oscs`--------------------------------
// Sinusoidal oscillator based on the state variable filter
// = undamped "modified-coupled-form" resonator
// = "magic circle" algorithm used in graphics.
//
// #### Usage
//
// ```
// oscs(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
oscs(f) = (*(-1) : sint(wn) : sintp(wn,impulse)) ~ _
with {
  wn = 2*ma.PI*f/ma.SR; // approximate
  // wn = 2*sin(PI*f/SR); // exact
  sint(x) = *(x) : + ~ _ ; // frequency-scaled integrator
  sintp(x,y) = *(x) : +(y): + ~ _; // same + state input
};

//-----------------`(os.)quadosc`--------------------
// Quadrature (cosine and sine) oscillator based on QuadOsc by Martin Vicanek.
//
// #### Usage
//
// ```
// quadosc(freq) : _,_
// ```
//
// where
//
// * `freq`: frequency in Hz
//
// #### Reference
// * <https://vicanek.de/articles/QuadOsc.pdf>
//------------------------------------------------------------
// Authors:
// Dario Sanfilippo <sanfilippo.dario@gmail.com>
// and Oleg Nesterov (JOS ed.)
quadosc(f) = tick ~ (_,_)
with {
    k1 = tan(f * ma.PI / ma.SR);
    k2 = 2 * k1 / (1 + k1 * k1);
    tick(u_0,v_0) = u_1,v_1
    with {
        tmp = u_0 - k1 * v_0;
        v_1 = v_0 + k2 * tmp;
        u_1 = tmp - k1 * v_1 : select2(1',1);
    };
};

//-----------------------------`(os.)sidebands`--------------------------------------
// Adds harmonics to quad oscillator.
//
// #### Usage
//
// ```
//    cos(x),sin(x) : sidebands(vs) : _,_
// ```
//
// Where:
//
// * `vs` : list of amplitudes
//
// #### Example test program
//
// ```
//    cos(x),sin(x) : sidebands((10,20,30))
// ```
//
// outputs:
//
// ```
//    10*cos(x) + 20*cos(2*x) + 30*cos(3*x),
//    10*sin(x) + 20*sin(2*x) + 30*sin(3*x);
// ```
//
// The following:
//
// ```
//    process = os.quadosc(F) : sidebands((10,20,30))
// ```
//
// is (modulo floating point issues) the same as:
//
// ```
//    c = os.quadosc : _,!;
//    s = os.quadosc : !,_;
//    process =
//        10*c(F) + 20*c(2*F) + 30*c(F),
//        10*s(F) + 20*s(2*F) + 30*s(F);
// ```
//
// but much more efficient.
//
// #### Implementation Notes
//
// This is based on the trivial trigonometric identities:
//
// ```
//    cos((n + 1) x) = 2 cos(x) cos(n x) - cos((n - 1) x)
//    sin((n + 1) x) = 2 cos(x) sin(n x) - sin((n - 1) x)
// ```
//
// Note that the calculation of the cosine/sine parts do not depend
// on each other, so if you only need the sine part you can do:
//
// ```
//    process = os.quadosc(F) : sidebands(vs) : !,_;
// ```
//
// and the compiler will discard the half of the calculations.
//-----------------------------------------------------------------------------
sidebands(vs, c0,s0)
	= c0*vn(0),s0*vn(0), 1,c0, 0,s0
	: seq(n, outputs(vs)-1, add(vn(n+1)))
	: _,_, !,!, !,!
with {
	// ba.take(n+1, vs)
	vn(n) = vs : route(outputs(vs),1, n+1,1);

	add(vn, co,so, cn_2,cn_1, sn_2,sn_1) =
		co+cn*vn, so+sn*vn, cn_1,cn, sn_1,sn
	with {
		cn = 2*c0*cn_1 - cn_2;
		sn = 2*c0*sn_1 - sn_2;
	};
};

//-----------------------------`(os.)sidebands_list`--------------------------------------
// Creates the list of complex harmonics from quad oscillator.
//
// Similar to `sidebands` but doesn't sum the harmonics, so it is more
// generic but less convenient for immediate usage.
//
// #### Usage
//
// ```
//    cos(x),sin(x) : sidebands_list(N) : si.bus(2*N)
// ```
//
// Where:
//
// * `N` : number of harmonics, compile time constant > 1
//
// #### Example test program
//
// ```
//    cos(x),sin(x) : sidebands_list(3)
// ```
//
// outputs:
//
// ```
//    cos(x),sin(x), cos(2*x),sin(2*x), cos(3*x),sin(3*x);
// ```
//
// The following:
//
// ```
//    process = os.quadosc(F) : sidebands_list(3)
// ```
//
// is (modulo floating point issues) the same as:
//
// ```
//    process = os.quadosc(F), os.quadosc(2*F), os.quadosc(3*F);
// ```
//
// but much more efficient.
//-----------------------------------------------------------------------------
sidebands_list(N, c0,s0)
	= c0,s0, 1,c0, 0,s0
	: seq(n, N-1, si.bus(2*(n+1)), add)
	: si.bus(2*N), !,!, !,!
with {
	add(cn_2,cn_1, sn_2,sn_1) =
		cn,sn, cn_1,cn, sn_1,sn
	with {
		cn = 2*c0*cn_1 - cn_2;
		sn = 2*c0*sn_1 - sn_2;
	};
};

//------------------------------`(os.)dsf`--------------------------------
// An environment with sine/cosine oscsillators with exponentially decaying
// harmonics based on direct summation formula.
//
// #### Usage
//
// ```
// dsf.xxx(f0, df, a, [n]) : _
// ```
//
// Where:
//
// * `f0`: base frequency
// * `df`: step frequency
// * `a`: decaying factor != 1
// * `n`: total number of harmonics (`osccN/oscsN` only)
//
// #### Variants
//
// - infinite number of harmonics, implies aliasing
// ```
// oscc(f0,df,a) : _;
// oscs(f0,df,a) : _;
// ```
//
// - n harmonics, f0, f0 + df, f0 + 2\*df, ..., f0 + (n-1)\*df
// ```
// osccN(f0,df,a,n) : _;
// oscsN(f0,df,a,n) : _;
// ```
//
// - finite number of harmonics, from f0 to Nyquist
// ```
// osccNq(f0,df,a) : _;
// oscsNq(f0,df,a) : _;
// ```
// 
// #### Example test program
//
// ```
// process = dsf.osccN(F0,DF,A,N),
//           dsf.oscsN(F0,DF,A,N);
// ```
// if `N` is an integer constant, the same (modulo fp issues) as:
// ```
// c = os.quadosc : _,!;
// s = os.quadosc : !,_;
// process = sum(k,N, A^k * c(F0 + k*DF)),
//           sum(k,N, A^k * s(F0 + k*DF));
// ```
// but much more efficient.
//
// #### Reference
//
// * <https://ccrma.stanford.edu/STANM/stanms/stanm5/stanm5.pdf>
//
declare dsf author "Oleg Nesterov";
dsf = environment {
	qo = quadosc;
	co = qo : _,!;
	so = qo : !,_;

	Nq(f0,df) = (ma.SR/2 - f0) / df : int;

	// infinite number of harmonics, implies aliasing
	oscc(f0,df,a) = (co(f0) - a*co(f0-df)) / (1 + a^2 - 2*a*co(df));
	oscs(f0,df,a) = (so(f0) - a*so(f0-df)) / (1 + a^2 - 2*a*co(df));

	// n harmonics, f0, f0 + df, f0 + 2*df, ..., f0 + (n-1)*df
	osccN(f0,df,a,n) = oscc(f0,df,a) - a^n * oscc(f0 + n*df, df,a);
	oscsN(f0,df,a,n) = oscs(f0,df,a) - a^n * oscs(f0 + n*df, df,a);

	// finite number of harmonics, from f0 to Nyquist
	osccNq(f0,df,a) = osccN(f0,df,a, Nq(f0,df));
	oscsNq(f0,df,a) = oscsN(f0,df,a, Nq(f0,df));
};

//================ Waveguide-Resonator-Based Oscillators ================
// Sinusoidal oscillator based on the waveguide resonator `wgr`.
//=======================================================================

//-----------------`(os.)oscwc`--------------------
// Sinusoidal oscillator based on the waveguide resonator `wgr`. Unit-amplitude
// cosine oscillator.
//
// #### Usage
//
// ```
// oscwc(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html>
//------------------------------------------------------------
oscwc(fr) = impulse : fi.wgr(fr,1) : _,!; // cosine (cheapest at 1 mpy/sample)

//-----------------`(os.)oscws`--------------------
// Sinusoidal oscillator based on the waveguide resonator `wgr`. Unit-amplitude
// sine oscillator.
//
// #### Usage
//
// ```
// oscws(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html>
//------------------------------------------------------------
oscws(fr) = impulse : fi.wgr(fr,1) : !,_; // sine (needs a 2nd scaling mpy)

//-----------------`(os.)oscq`--------------------
// Sinusoidal oscillator based on the waveguide resonator `wgr`.
// Unit-amplitude cosine and sine (quadrature) oscillator.
//
// #### Usage
//
// ```
// oscq(freq) : _,_
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html>
//------------------------------------------------------------
oscq(fr) = impulse : fi.wgr(fr,1);       // phase quadrature outputs

//-----------------`(os.)oscw`--------------------
// Sinusoidal oscillator based on the waveguide resonator `wgr`.
// Unit-amplitude cosine oscillator (default).
//
// #### Usage
//
// ```
// oscw(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html>
//------------------------------------------------------------
oscw = oscwc;

// end jos section
//########################################################################################
/************************************************************************
FAUST library file, further contributions section

All contributions below should indicate both the contributor and terms
of license.  If no such indication is found, "git blame" will say who
last edited each line, and that person can be emailed to inquire about
license disposition, if their license choice is not already indicated
elsewhere among the libraries.  It is expected that all software will be
released under LGPL, STK-4.3, MIT, BSD, or a similar FOSS license.
************************************************************************/

//===================== Casio CZ Oscillators ==========================
// Oscillators that mimic some of the Casio CZ oscillators.
//
// There are two sets:
//
// * a set with an index parameter
//
// * a set with a res parameter
//
// The "index oscillators" outputs a sine wave at index=0 and gets brighter with a higher index.
// There are two versions of the "index oscillators":
//
// * with P appended to the name: is phase aligned with `fund:sin`
//
// * without P appended to the name: has the phase of the original CZ oscillators
//
// The "res oscillators" have a resonant frequency.
// "res" is the frequency of resonance as a factor of the fundamental pitch.
//
// For the `fund` waveform, use a low-frequency oscillator without anti-aliasing such as `os.lf_saw`.
//=====================================================================

//----------`(os.)CZsaw`----------
// Oscillator that mimics the Casio CZ saw oscillator.
// `CZsaw` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsaw(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 to 1. 0 = sine-wave, 1 = saw-wave
//------------------------------------------------------------
declare CZsaw author "Bart Brouns";
declare CZsaw licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsaw(fund, index) = CZ.sawChooseP(fund, index, 0);

//----------`(os.)CZsawP`----------
// Oscillator that mimics the Casio CZ saw oscillator,
// with it's phase aligned to `fund:sin`.
// `CZsawP` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsawP(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 to 1. 0 = sine-wave, 1 = saw-wave
//------------------------------------------------------------
declare CZsawP author "Bart Brouns";
declare CZsawP licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsawP(fund, index) = CZ.sawChooseP(fund, index, 1);

//----------`(os.)CZsquare`----------
// Oscillator that mimics the Casio CZ square oscillator
// `CZsquare` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsquare(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 to 1. 0 = sine-wave, 1 = square-wave
//------------------------------------------------------------
declare CZsquare author "Bart Brouns";
declare CZsquare licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsquare(fund, index) = CZ.squareChooseP(fund, index, 0);

//----------`(os.)CZsquareP`----------
// Oscillator that mimics the Casio CZ square oscillator,
// with it's phase aligned to `fund:sin`.
// `CZsquareP` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsquareP(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 to 1. 0 = sine-wave, 1 = square-wave
//------------------------------------------------------------
declare CZsquareP author "Bart Brouns";
declare CZsquareP licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsquareP(fund, index) = CZ.squareChooseP(fund, index, 1);

//----------`(os.)CZpulse`----------
// Oscillator that mimics the Casio CZ pulse oscillator.
// `CZpulse` is a standard Faust function.
//
// #### Usage
//
// ```
// CZpulse(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is closer to a pulse
//------------------------------------------------------------
declare CZpulse author "Bart Brouns";
declare CZpulse licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZpulse(fund, index) = CZ.pulseChooseP(fund, index, 0);

//----------`(os.)CZpulseP`----------
// Oscillator that mimics the Casio CZ pulse oscillator,
// with it's phase aligned to `fund:sin`.
// `CZpulseP` is a standard Faust function.
//
// #### Usage
//
// ```
// CZpulseP(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is closer to a pulse
//------------------------------------------------------------
declare CZpulseP author "Bart Brouns";
declare CZpulseP licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZpulseP(fund, index) = CZ.pulseChooseP(fund, index, 1);

//----------`(os.)CZsinePulse`----------
// Oscillator that mimics the Casio CZ sine/pulse oscillator.
// `CZsinePulse` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsinePulse(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is a sine minus a pulse
//------------------------------------------------------------
declare CZsinePulse author "Bart Brouns";
declare CZsinePulse licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsinePulse(fund, index) = CZ.sinePulseChooseP(fund, index, 0);

//----------`(os.)CZsinePulseP`----------
// Oscillator that mimics the Casio CZ sine/pulse oscillator,
// with it's phase aligned to `fund:sin`.
// `CZsinePulseP` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsinePulseP(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is a sine minus a pulse
//------------------------------------------------------------
declare CZsinePulseP author "Bart Brouns";
declare CZsinePulseP licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsinePulseP(fund, index) = CZ.sinePulseChooseP(fund, index, 1);

//----------`(os.)CZhalfSine`----------
// Oscillator that mimics the Casio CZ half sine oscillator.
// `CZhalfSine` is a standard Faust function.
//
// #### Usage
//
// ```
// CZhalfSine(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is somewhere between a saw and a square
//------------------------------------------------------------
declare CZhalfSine author "Bart Brouns";
declare CZhalfSine licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZhalfSine(fund, index) = CZ.halfSineChooseP(fund, index, 0);

//----------`(os.)CZhalfSineP`----------
// Oscillator that mimics the Casio CZ half sine oscillator,
// with it's phase aligned to `fund:sin`.
// `CZhalfSineP` is a standard Faust function.
//
// #### Usage
//
// ```
// CZhalfSineP(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is somewhere between a saw and a square
//------------------------------------------------------------
declare CZhalfSineP author "Bart Brouns";
declare CZhalfSineP licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZhalfSineP(fund, index) = CZ.halfSineChooseP(fund, index, 1);

//----------`(os.)CZresSaw`----------
// Oscillator that mimics the Casio CZ resonant sawtooth oscillator.
// `CZresSaw` is a standard Faust function.
//
// #### Usage
//
// ```
// CZresSaw(fund,res) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `res`: the frequency of resonance as a factor of the fundamental pitch.
//------------------------------------------------------------
declare CZresSaw author "Bart Brouns";
declare CZresSaw licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZresSaw(fund,res) = CZ.resSaw(fund,res);

//----------`(os.)CZresTriangle`----------
// Oscillator that mimics the Casio CZ resonant triangle oscillator.
// `CZresTriangle` is a standard Faust function.
//
// #### Usage
//
// ```
// CZresTriangle(fund,res) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `res`: the frequency of resonance as a factor of the fundamental pitch.
//------------------------------------------------------------
declare CZresTriangle author "Bart Brouns";
declare CZresTriangle licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZresTriangle(fund,res) = CZ.resTriangle(fund,res);

//----------`(os.)CZresTrap`----------
// Oscillator that mimics the Casio CZ resonant trapeze oscillator
// `CZresTrap` is a standard Faust function.
//
// #### Usage
//
// ```
// CZresTrap(fund,res) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `res`: the frequency of resonance as a factor of the fundamental pitch.
//------------------------------------------------------------
declare CZresTrap author "Bart Brouns";
declare CZresTrap licence "STK-4.3";

// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZresTrap(fund, res) = CZ.resTrap(fund, res);

CZ = environment {

    saw(fund, index) = sawChooseP(fund, index, 0);
    sawP(fund, index) = sawChooseP(fund, index, 1);
    sawChooseP(fund, index, p) =
      (((FUND(fund,align,p)*((.5-INDEX)/INDEX)),(-1*FUND(fund,align,p)+1)*((.5-INDEX)/(1-INDEX))):min+FUND(fund,align,p))*2*ma.PI:cos
    with {
      INDEX = (.5-(index*.5)):max(0.01):min(0.5);
      align = si.interpolate(index, 0.75, 0.5);
    };

    square(fund, index) = squareChooseP(fund, index, 0);
    squareP(fund, index) = squareChooseP(fund, index, 1);
    squareChooseP(fund, index, p) = (FUND(fund,align,p)>=0.5), (ma.decimal((FUND(fund,align,p)*2)+1)<:_-min(_,(-1*_+1)*((INDEX)/(1-INDEX)))) :+ *ma.PI:cos
    with {
      INDEX = (index:pow(0.25)):max(0):min(1);
      align = si.interpolate(INDEX, -0.25, 0);
    };

    pulse(fund, index) = pulseChooseP(fund, index, 0);
    pulseP(fund, index) = pulseChooseP(fund, index, 1);
    pulseChooseP(fund, index, p) = ((FUND(fund,align,p)-min(FUND(fund,align,p),((-1*FUND(fund,align,p)+1)*(INDEX/(1-INDEX)))))*2*ma.PI):cos
    with {
      INDEX = index:min(0.99):max(0);
      align = si.interpolate(index, -0.25, 0.0);
    };

    sinePulse(fund, index) = sinePulseChooseP(fund, index, 0);
    sinePulseP(fund, index) = sinePulseChooseP(fund, index, 1);
    sinePulseChooseP(fund, index, p) = (min(FUND(fund,align,p)*((0.5-INDEX)/INDEX),(-1*FUND(fund,align,p)+1)*((.5-INDEX)/(1-INDEX)))+FUND(fund,align,p))*4*ma.PI:cos
    with {
      INDEX = ((index*-0.49)+0.5);
      align = si.interpolate(index, -0.125, -0.25);
    };

    halfSine(fund, index) = halfSineChooseP(fund, index, 0);
    halfSineP(fund, index) = halfSineChooseP(fund, index, 1);
    halfSineChooseP(fund, index, p) = (select2(FUND(fund,align,p)<.5, .5*(FUND(fund,align,p)-.5)/INDEX+.5, FUND(fund,align,p)):min(1))*2*ma.PI:cos
    with {
      INDEX = (.5-(index*0.5)):min(.5):max(.01);
      align = si.interpolate(index:min(0.975), -0.25, -0.5);
    };

    FUND =
      case {
        (fund,align,0) => fund;
        (fund,align,1) => (fund+align) : ma.frac; // align phase with fund
      };
    resSaw(fund,res) = (((-1*(1-fund))*((cos((ma.decimal((max(1,res)*fund)+1))*2*ma.PI)*-.5)+.5))*2)+1;
    resTriangle(fund,res) = select2(fund<.5, 2-(fund*2), fund*2)*INDEX*2-1
    with {
      INDEX = ((fund*(res:max(1)))+1:ma.decimal)*2*ma.PI:cos*.5+.5;
    };
    resTrap(fund, res) = (((1-fund)*2):min(1)*sin(ma.decimal(fund*(res:max(1)))*2*ma.PI));
};

//===============================PolyBLEP-Based Oscillators=================================

//----------`(os.)polyblep`----------
// PolyBLEP residual function, used for smoothing steps in the audio signal.
//
// #### Usage
//
// ```
// polyblep(Q,phase) : _
// ```
//
// Where:
//
// * `Q`: smoothing factor between 0 and 0.5. Determines how far from the ends of the phase interval the quadratic function is used.
// * `phase`: normalised phase (between 0 and 1)
//------------------------------------------------------------
declare polyblep author "Jacek Wieczorek";

polyblep(Q, phase) = (0, L(phase / Q), R((phase - 1) / Q)) : select3(sel)
with {
    sel = (phase < Q) + 2*(phase > 1 - Q);
    L(x) = 2*x - x*x - 1; // Used near the left end of the interval
    R(x) = 2*x + x*x + 1; // Used near the right end of the interval
};

//----------`(os.)polyblep_saw`----------
// Sawtooth oscillator with suppressed aliasing (using `polyblep`).
//
// #### Usage
//
// ```
// polyblep_saw(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
declare polyblep_saw author "Jacek Wieczorek";

polyblep_saw(freq) = naive - polyblep(Q , phase)
with {
    phase = phasor(1, freq);
    naive = 2 * phase - 1;
    Q = freq / ma.SR;
};

//----------`(os.)polyblep_square`----------
// Square wave oscillator with suppressed aliasing (using `polyblep`).
//
// #### Usage
//
// ```
// polyblep_square(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
declare polyblep_square author "Jacek Wieczorek";

polyblep_square(freq) = naive - polyblep(Q, phase) + polyblep(Q, ma.modulo(phase + 0.5, 1))
with {
    phase = phasor(1, freq);
    naive = 2 * (phase * 2 : int) - 1;
    Q = freq / ma.SR;
};

//----------`(os.)polyblep_triangle`----------
// Triangle wave oscillator with suppressed aliasing (using `polyblep`).
//
// #### Usage
//
// ```
// polyblep_triangle(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
declare polyblep_triangle author "Jacek Wieczorek";

polyblep_triangle(freq) = polyblep_square(freq) : fi.pole(0.999) : *(4 * freq / ma.SR);

// end further contributions section
//########################################################################################