declare name "Faust Oscillator Library";
declare author "Julius O. Smith (jos at ccrma.stanford.edu)";
declare copyright "Julius O. Smith III";
declare version "1.10";
declare license "STK-4.3"; // Synthesis Tool Kit 4.3 (MIT style license)

import("music.lib");  // SR, ...
import("filter.lib"); // wgr, nlf2, tf2

//===================== Virtual Analog Oscillators ========================

//------------------------ Impulse Train: imptrain ------------------------
imptrain(freq) = sawpos(freq)<:-(mem)<0;

//--- Pulse-Train and Square-Wave Oscillators: pulsetrainpos, squarewave[pos]
// In all cases, the first pulse jumps to 1 at time 0.

// Basic unit-amplitude nonnegative pulse train with duty cycle between 0 and 1:
pulsetrainpos(freq,duty) = float(sawpos(freq) <= duty);

// Positive square wave = pulse train with 50% duty cycle:
squarewavepos(freq) = pulsetrainpos(freq,0.5);

// Unit amplitude square wave = zero-mean pulse train with 50% duty cycle:
squarewave(freq) = 2*squarewavepos(freq) - 1;

//---------- Sawtooth: rawsaw, sawpos, saw1, saw2, sawtooth -------------

// Sawtooth waveform oscillators for virtual analog synthesis et al.
// The 'simple' versions (rawsaw, sawpos, 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, 

// --- rawsaw ---
// simple sawtooth waveform oscillator between 0 and period in samples:
rawsaw(periodsamps) = (_,periodsamps : fmod) ~ +(1.0);

// --- sawpos ---
// simple sawtooth waveform oscillator between 0 and 1
sawpos(freq) = rawsaw(periodsamps) / periodsamps
with {
  periodsamps = float(SR)/freq; // period in samples (not nec. integer)
};

// --- saw1 ---
// simple sawtooth waveform oscillator between -1 and 1
saw1(freq) = 2.0 * sawpos(freq) - 1.0; // zero-mean in [-1,1)

// --- saw2 ---
// Differentiated Parabolic Wave sawtooth (less aliasing)
// Reference: Valimaki, IEEE Signal Processing Letters, March 2005
saw2(freq) = saw1(freq) <: * <: -(mem) : *(0.25'*SR/freq);

// --- sawtooth ---
sawtooth = saw2; // default choice

//-------------------------- sawtooth_demo ---------------------------
// USAGE:  sawtooth_demo : _

sawtooth_demo = signal with { 
  osc_group(x) = vgroup("[0] SAWTOOTH OSCILLATOR 
    [tooltip: See Faust's oscillator.lib for documentation and references]",x);
  knob_group(x)  = osc_group(hgroup("[1]", x));
  ampdb  = knob_group(vslider("[1] Amplitude [unit:dB] [style:knob]
    [tooltip: Sawtooth waveform amplitude]",
    -20,-120,10,0.1));
  amp = ampdb : smooth(0.999) : db2linear;
  freq = knob_group(
  vslider("[2] Frequency [unit:PK] [style:knob]
    [tooltip: Sawtooth frequency as a Piano Key (PK) number (A440 = key 49)]",
    49,1,88,0.01) : pianokey2hz);
  pianokey2hz(x) = 440.0*pow(2.0, (x-49.0)/12); // piano key 49 = A440 (also defined in effect.lib)
  detune1 = 1 + 0.01 * knob_group(
    vslider("[3] Detuning 1 [unit:%%] [style:knob]
      [tooltip: Percentange frequency-shift up or down for second oscillator]",
      -0.1,-10,10,0.01));
  detune2 = 1 + 0.01 * knob_group(
    vslider("[4] Detuning 2 [unit:%%] [style:knob]
[tooltip: Percentange frequency-shift up or down for third detuned oscillator]",
    +0.1,-10,10,0.01));
  portamento = knob_group(
    vslider("[5] Portamento [unit:sec] [style:knob]
      [tooltip: Portamento (frequency-glide) time-constant in seconds]",
      0.1,0.01,1,0.001));
  sfreq = freq : smooth(tau2pole(portamento));
  tone = (amp/3) * 
    (sawtooth(sfreq) + sawtooth(sfreq*detune1) + sawtooth(sfreq*detune2));
  signal = amp * select2(ei, select2(ss, tone, pink_noise), _);
  checkbox_group(x) = knob_group(vgroup("[6] Alternate Signals",x));
  ss = checkbox_group(checkbox("[0] 
[tooltip: Pink Noise (or 1/f noise) is Constant-Q Noise, meaning that it has the same total power in every octave] Pink Noise Instead (uses only Amplitude control on the left)"));
  ei = checkbox_group(checkbox(
      "[1] External Input Instead (overrides Sawtooth/Noise selection above)"));
};

// --- Correction-filtered versions of saw2: saw2f2, saw2f4 ----
saw2f2 = saw2 : cf2 with {
  cf2 = tf2(1.155704605878911, 0.745184288225518,0.040305967265900,
	0.823765146386639, 0.117420665547108);
};

saw2f4 = saw2 : cf4 with {
  cf4 = iir((1.155727435125014, 2.285861038554662,
	1.430915027294021, 0.290713280893317, 0.008306401748854),
	(2.156834679164532, 1.559532244409321, 0.423036498118354,
	0.032080681130972));
};

// --- sawN, saw3,saw4,saw5,saw6 ---
// Differentiated Polynomial Wave (DPW) sawtooth (progressively less aliasing)
// Reference: 
// "Alias-Suppressed Oscillators based on Differentiated Polynomial Waveforms",
// Vesa Valimaki, Juhan Nam, Julius Smith, and Jonathan Abel, 
// IEEE Tr. Acoustics, Speech, and Language Processing (IEEE-ASLP),
// Vol. 18, no. 5, May 2010.

sawN(N,freq) = saw1 : poly(N) : D(N-1) : gate(N-1)
with {
  p0n = SR/freq;
  sawpos = (_,1:fmod) ~ +(1/p0n); // sawtooth waveform in [0,1)
  saw1 = 2*sawpos - 1; // zero average mean, unit max amp
  poly(2,x) =  x*x;
  poly(3,x) =  x*x*x - x;
  poly(4,x) =  poly(2,x)*(poly(2,x) - 2);
  poly(5,x) =  pow(x,5) - pow(x,3)*10/3 + x*7/3;
  poly(6,x) =  pow(x,6) - 5*pow(x,4) + 7*poly(2,x);
  diff1(x) =  (x - x')/(2/p0n);
  diff(N) = seq(n,N,diff1); // N diffs in series
  D(1) = diff1/2;
  D(2) = diff(2)/6;
  D(3) = diff(3)/24;
  D(4) = diff(4)/120;
  D(5) = diff(5)/720;
  gatedelay(n,d,x)    = x@(int(d)&(n-1)); // from music.lib
  gate(N) = * (1 : gatedelay(8,N)); // delayed step for blanking startup glitch
};
saw3 = sawN(3); saw4 = sawN(4); saw5 = sawN(5); saw6 = sawN(6);

//----------------------- Filter-Based Oscillators ------------------------

// Quick Guide (more complete documentation forthcoming):
//
//    USAGE: osc[b|r|rs|rc|s|w](f), where f = frequency in Hz.
//
//  oscb: one-multiply, two-adds, amplitude varies with frequency, avoid dc
//  oscr: four-multipies, two-adds, amplitude unchanging with frequency, 
//        dc ok, amp slowly drifts,
//        sine and cosine outputs available (exact phase quadrature)
// oscrs: sine output of oscr
// oscrc: cosine output of oscr
//  oscs: two-multiplies, two-adds, amplitude varies slightly with frequency, 
//        dc ok, no amp drift, likely optimizable to be the fastest no-drift case
//  oscw: one/two-multiply, three-adds, amplitude steady with frequency, no amp drift,
//        sine and cosine outputs available (exact phase quadrature),
//        numerical difficulty below 10 Hz,
//        likely optimizable to be best (above 10 Hz) for custom silicon
//        (one multiply when frequency is constant, two otherwise).

impulse = 1-1'; // used to start filter-based oscillators

//-------------------------- oscb --------------------------------
// Sinusoidal oscillator based on the biquad
//
oscb(f) = impulse : tf2(1,0,0,a1,1)
with {
  a1 = -2*cos(2*PI*f/SR);
};

//-------------------------- oscr --------------------------------
// Sinusoidal oscillator based on 2D vector rotation,
//  = undamped "coupled-form" resonator
//  = lossless 2nd-order normalized ladder filter
//
// Reference:
// https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html
//
oscrq(f) =  impulse : nlf2(f,1); // sine and cosine outputs
oscrs(f) =  impulse : nlf2(f,1) : _,!; // sine
oscrc(f) =  impulse : nlf2(f,1) : !,_; // cosine
oscr = oscrs; // default = sine case

//-------------------------- oscs --------------------------------
// Sinusoidal oscillator based on the state variable filter
// = undamped "modified-coupled-form" resonator
//
oscs(f) =  (*(0-1) : sint(wn) : sintp(wn,impulse)) ~ _
with {
  wn = 2*PI*f/SR; // approximate
  // wn = 2*sin(PI*f/SR); // exact
  sub(x,y) = y-x;
  sint(x) = *(x) : + ~ _ ; // frequency-scaled integrator
  sintp(x,y) = *(x) : +(y): + ~ _ ; // same + state input
};

//----------------- oscw, oscwq, oscwc, oscws --------------------
// Sinusoidal oscillator based on the waveguide resonator wgr
//
// oscwc - unit-amplitude cosine oscillator
// oscws - unit-amplitude sine oscillator
// oscq  - unit-amplitude cosine and sine (quadrature) oscillator
// oscw  - default = oscwc for maximum speed
//
// Reference:
// https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html
//
oscwc(fr) = impulse : wgr(fr,1) : _,!; // cosine (cheapest at 1 mpy/sample)
oscws(fr) = impulse : wgr(fr,1) : !,_; // sine (needs a 2nd scaling mpy)
oscq(fr)  = impulse : wgr(fr,1);       // phase quadrature outputs
oscw = oscwc;

//-------------------------- oscrs_demo ---------------------------

oscrs_demo = signal with { 
  osc_group(x) = vgroup("[0] SINE WAVE OSCILLATOR oscrs 
    [tooltip: Sine oscillator based on 2D vector rotation]",x);
  knob_group(x)  = osc_group(hgroup("[1]", x));
//  ampdb  = knob_group(vslider("[1] Amplitude [unit:dB] [style:knob] 
  ampdb  = knob_group(hslider("[1] Amplitude [unit:dB]
    [tooltip: Sawtooth waveform amplitude]",
    -20,-120,10,0.1));
  amp = ampdb : smooth(0.999) : db2linear;
  freq = knob_group(
//  vslider("[2] Frequency [unit:PK] [style:knob] 
  hslider("[2] Frequency [unit:PK]
     [tooltip: Sine wave frequency as a Piano Key (PK) number (A440 = 49 PK)]",
     49,1,88,0.01) : pianokey2hz);
  pianokey2hz(x) = 440.0*pow(2.0, (x-49.0)/12); // (also defined in effect.lib)
  portamento = knob_group(
//    vslider("[3] Portamento [unit:sec] [style:knob] 
    hslider("[3] Portamento [unit:sec]
      [tooltip: Portamento (frequency-glide) time-constant in seconds]",
      0.1,0,1,0.001));
  sfreq = freq : smooth(tau2pole(portamento));
  signal = amp * oscrs(sfreq);
};

oscr_demo = oscrs_demo; // synonym

//--------------------------- pink_noise --------------------------
// Pink noise (1/f noise) generator (third-order approximation)
//
// USAGE: pink_noise : _;
//
// Reference:
//  https://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html
//

pink_noise = noise : 
   iir((0.049922035, -0.095993537, 0.050612699, -0.004408786),
                    (-2.494956002, 2.017265875, -0.522189400));