desc:Saike Tanh Saturation with anti aliasing
tags: saturation distortion anti-aliased
version: 1.16
author: Joep Vanlier, Erich M. Burg
changelog: 
  + Added oversampling
provides: upsampler/*
license: MIT

Uses technique from: Parker et al, "REDUCING THE ALIASING OF NONLINEAR WAVESHAPING USING CONTINUOUS-TIME CONVOLUTION",
Proceedings of the 19th International Conference on Digital Audio Effects (DAFx-16), Brno, Czech Republic, September 5–9, 2016
Special thanks to Erich M. Burg for coming up with the odd function mirror trick to improve numerical stability and
parameter smoothing.

import tanh_upsampler.jsfx-inc

in_pin:left input
in_pin:right input
out_pin:left output
out_pin:right output

slider1:0<-6,24,1>Gain (dB)
slider2:0<-18,0,1>Ceiling (dB)
slider3:1<0,2,1{No,Constant,Linear (BETA)}>Continuous antialias mode
slider4:0<0,1,1>Fix DC?
slider5:smooth_time=50<0,100,.001>Inertia (ms)
slider6:shelf=1<0,1,1>Correct for HF loss
slider7:func=0<0,3,1{tanh,clip,sin}>Shaping function
slider8:Oversampling=1<1,4,1>Oversampling

@init
  function generate_shelf()
  local(alpha, beta, twodivpisq, f, ifcsq, theta, b, G)
  global()
  instance(a1, a2, b0, b1)
  (
    G = 2.20;
    f = 2.5;
    theta = 1.0 - cos($pi * f);
    twodivpisq = 2.0 / $pi * $pi;
    
    ifcsq = 1.0/(f*f);
    alpha = twodivpisq * ( ifcsq + ifcsq / G ) - theta;
    beta = twodivpisq * (ifcsq + G * ifcsq ) - theta;
    
    a1 = - alpha / ( 1.0 + alpha + sqrt(1.0 + 2.0*alpha) );
    b = - beta / ( 1.0 + beta + sqrt(1.0 + 2.0*beta) );
    b0 = (1.0 + a1) / (1.0 + b);
    b1 = b * b0;
  );
  
  function generate_shelf_alternate()
  local(alpha, beta, twodivpisq, f, ifcsq, theta, b, G)
  global()
  instance(a1, a2, b0, b1)
  (
    G = 2.20;
    f = 2.5;
    theta = 1.0 - cos($pi * f);
    twodivpisq = 2.0 / $pi * $pi;
    
    ifcsq = 1.0/(f*f);
    alpha = twodivpisq * ( ifcsq + ifcsq / G ) - theta;
    beta = twodivpisq * (ifcsq + G * ifcsq ) - theta;
    
    alpha = twodivpisq * ( 1 + ifcsq / G ) - .5;
    beta = twodivpisq * ( 1 + G * ifcsq ) - .5;
        
    a1 = - alpha / ( 1.0 + alpha + sqrt(1.0 + 2.0*alpha) );
    b = - beta / ( 1.0 + beta + sqrt(1.0 + 2.0*beta) );
    b0 = (1.0 + a1) / (1.0 + b);
    b1 = b * b0;
  );  

 function smooth_parameter(target)
  instance(s, coeff)
  global(maxo, mino)
  local(y, v)
  (
    v = coeff*(target - s);
    y = v + s;
    s = y + v;
    y
  );
  
  function initialize_smoother(cutoff)
  instance(coeff)
  global(srate)
  local(g)
  ( 
    g = tan($pi*cutoff/srate);
    coeff = g/(1+g);  
  );

  function interpolator_init(slider_idx)
  instance(next_val)
  local()
  global()
  (
    next_val = slider(slider_idx);
  );

  function interpolator_block(slider_idx)
  instance(delta, next_changepoint_t, next_val)
  local(next_changepoint_y)
  global(samplesblock)
  (
    next_changepoint_t = slider_next_chg(slider_idx, next_changepoint_y);
    next_changepoint_t > 0 ? 
    (
      next_val = slider(slider_idx);
    ) : (
      next_changepoint_y = slider(slider_idx);
      next_changepoint_t = samplesblock;
    );
    
    delta = (next_changepoint_y - next_val) / next_changepoint_t;
  );
  
  function interpolate(slider_idx)
  instance(delta, next_changepoint_t, next_val)
  local(current_value, next_changepoint_y)
  global(current_sample)
  (
    current_value = next_val;
    current_sample == next_changepoint_t ? 
    (
      delta = 0;
      next_changepoint_t = slider_next_chg(slider_idx, next_changepoint_y);
      delta = next_changepoint_t > current_sample ? (next_changepoint_y - current_value) / (next_changepoint_t-current_sample) : 0;
    );
    
    next_val = current_value + delta;
    current_value
  );
  
  function initialize_interpolators()
  (
    gain_interpolator.interpolator_init(1);
    ceiling_interpolator.interpolator_init(2);
    new_cutoff = 1000.0/(smooth_time);
    (new_cutoff != cutoff) ? (
      cutoff = new_cutoff;
      gain_smoother.initialize_smoother(cutoff);
      ceiling_smoother.initialize_smoother(cutoff);
      gain_smoother.s = slider1;
      ceiling_smoother.s = slider2;
    );
    l_shelf.generate_shelf();
    r_shelf.generate_shelf();
  );
  
  initialize_interpolators();

@slider
  initialize_interpolators();
  
  (slider3 == 2) ? (
    target_pdc = 1 + getFIRdelay(overSampling);
  ) : (slider3 == 1) ? (
    target_pdc = 1 + getFIRdelay(overSampling);
  ) : (slider3 == 0) ? (
    target_pdc = getFIRdelay(overSampling);
  );

@block
  target_pdc != current_pdc ? (
    pdc_bot_ch = 0;
    pdc_top_ch = 2;
    pdc_delay = target_pdc;
    current_pdc = target_pdc;
  );

  current_sample = 0;
  gain_interpolator.interpolator_block(1);
  ceiling_interpolator.interpolator_block(2);

@sample
  /* Sample accurate interpolation */
  current_gain = gain_interpolator.interpolate(1);
  current_ceiling = ceiling_interpolator.interpolate(2);
    
  smooth_time ? (
    current_gain = gain_smoother.smooth_parameter(current_gain);
    current_ceiling = ceiling_smoother.smooth_parameter(current_ceiling);
  );
  
  // exp((log(10)/20) * X) is equivalent but faster than 10^(X/20), which is why we use this form
  log10d20_conversion  = .11512925464970228420089957273422; //log(10)/20;
  preamp      = exp(log10d20_conversion * current_gain);
  ceiling     = exp(-log10d20_conversion * current_ceiling);
  inv_ceiling = 1.0 / ceiling;

  function F0_sin(x)
  (
    x > 1.0 ? x + 2.0/$pi - 1.0 : 
    x < -1.0 ? -x + 2.0/$pi - 1.0 :
    -2.0*(cos(0.5*$pi*x)-1.0)/$pi
  );

  function shape_sin(x)
  local()
  global()
  (
    x > 1.0 ? 1.0 : 
    x < -1.0 ? -1.0 :
    sin(.5*$pi*x)
  );

  function F0_clip(x)
  (
    x > 1.0 ? x - 0.5 : 
    x < -1.0 ? -x - 0.5 :
    .5*x*x
  );

  function clip(x)
  local()
  global()
  (
    x > 1.0 ? 1.0 : 
    x < -1.0 ? -1.0 :
    x
  );

  function F0_tanh(x)
  local()
  global()
  (
    x - log(2/(1 + exp(-2*x)))
  );
  
  function tanh(x)
  local()
  global()
  instance()
  (
    2/(1+exp(-2*x)) - 1
  );

  function Li2(x)
  local(A, ALFA, B0, B1, B2, H, S, T, Y, Q, HF, PI3, PI6, PI12)
  global()
  (
    HF = 0.5;
    PI3 = 3.2898681336964528729448303332921; // pi*pi/3
    PI6 = 1.644934066848226436472415166646; // pi*pi / 6
    PI12 = 0.82246703342411321823620758332301; // pi*pi / 12
    
    (x==1) ? (
      H = PI6;
    ) : (x == -1) ? (
      H = -PI12;
    ) : (
      T = -x;
    );

    A = (T <= -2) ? (
      Y = -1 / (1 + T);
      S = 1;
      A = log(-T);
      Q = log(1 + 1/T);
      -PI3 + HF * (A*A - Q*Q)
    ) : (T < -1) ? (
      Y = -1 - T;
      S = -1;
      A = log(-T);
      Q = log(1 + 1/T);  
      -PI6 + A * (A + Q)
    ) : (T <= -0.5) ? (
      Y = -(1 + T) / T;
      S = 1;
      A = log(-T);
      Q = log(1 + T);
      -PI6 + A * (-HF * A + Q)
    ) : (T < 0) ? (
      Y = -T / (1 + T);
      S = -1;
      Q = log(1 + T);
      A = HF * Q*Q
    ) : (T <= 1) ? (
      Y = T;
      S = 1;
      0
    ) : (
      Y = 1 / T;
      S = -1;
      Q = log(T);
      PI6 + HF * Q*Q
    );

    H = Y + Y - 1;
    ALFA = H + H;
    B1 = 0;
    B2 = 0;
    B0 = 0.00000000000000002 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = -0.00000000000000014 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = 0.00000000000000093 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = -0.00000000000000610 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = 0.00000000000004042 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = -0.00000000000027007 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = 0.00000000000182256 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = -0.00000000001244332 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = 0.00000000008612098 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = -0.00000000060578480 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = 0.00000000434545063 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = -0.00000003193341274 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = 0.00000024195180854 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = -0.00000190784959387 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = 0.00001588415541880 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = -0.00014304184442340 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = 0.00145751084062268 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = -0.01858843665014592 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = 0.40975987533077105 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    B0 = 0.42996693560813697 + ALFA * B1 - B2; B2 = B1; B1 = B0;
    
    H = -(S * (B0 - H * B2) + A)
  );
  
  function F1(x)
  local(em2x)
  global()
  (
    em2x = exp(-2*x);
    .5 * (x * (x + 2 * log(em2x + 1)) - Li2(-em2x))
  );
  
  function antialiased_tanh_linear(xn)
  local(eps, absxn, hpi12)
  global()
  instance(antialias, F0_xnm1, xnm1, xnm2, F0_xnm1, F0_xnm2, F1_xnm1, F1_xnm2,
  F0_xn, F1_xn, diff1, diff2, term1, term2, idiff)
  (
    absxn     = abs(xn);
    F0_xn     = F0_tanh(absxn);
    
    hpi12     = 0.4112335167120566;  // .5 * pi*pi / 12
    F1_xn     = (F1(absxn) - hpi12)*sign(xn) + hpi12;
 
    diff1     = ( xn - xnm1 );
    diff2     = ( xnm2 - xnm1 );
    eps       = .2;
      
    term1     = (abs(diff1) > eps) ? (
      idiff = 1.0 / (diff1*diff1);
      ( xn * ( F0_xn - F0_xnm1 ) - (F1_xn - F1_xnm1) ) * idiff
    ) : (
      .5 * tanh((xn + 2*xnm1)*.33333333333333333333333333333)
    );
    
    term2     = (abs(diff2) > eps) ? (
      idiff = 1.0 / (diff2*diff2);
      ( xnm2 * ( F0_xnm2 - F0_xnm1 ) - (F1_xnm2 - F1_xnm1) ) * idiff
    ) : (
      .5 * tanh((xnm2 + 2*xnm1)*.33333333333333333333333333333)
    );

    F1_xnm2   = F1_xnm1;
    F1_xnm1   = F1_xn;

    F0_xnm2   = F0_xnm1;
    F0_xnm1   = F0_xn;
    
    xnm2      = xnm1;
    xnm1      = xn;
  
    term1 + term2
  );
  
  function antialiased_tanh_rect(x)
  local(eps, F0_xn, diff)
  global()
  instance(antialias, F0_xnm1, xnm1)
  (
    (
      F0_xn     = F0_tanh(abs(x));
      diff      = ( x - xnm1 );
      eps       = 0.0000000001;
      antialias = (abs(diff) > eps) ? ( F0_xn - F0_xnm1 ) / diff : tanh(.5*(x+xnm1));
    );
    
    F0_xnm1   = F0_xn;
    xnm1      = x;

    antialias
  );
  
  function antialiased_clip_rect(x)
  local(eps, F0_xn)
  global()
  instance(antialias, F0_xnm1, xnm1,diff)
  (
    (
      F0_xn     = F0_clip(abs(x));
      diff      = ( x - xnm1 );
      eps       = 0.0000000001;
      antialias = (abs(diff) > eps) ? ( F0_xn - F0_xnm1 ) / diff : clip(.5*(x+xnm1));
    );
    
    F0_xnm1   = F0_xn;
    xnm1      = x;
  
    antialias
  );
  
  function antialiased_sin_rect(x)
  local(eps, F0_xn)
  global()
  instance(antialias, F0_xnm1, xnm1,diff)
  (
    (
      F0_xn     = F0_sin(x);
      diff      = ( x - xnm1 );
      eps       = 0.0000000001;
      antialias = (abs(diff) > eps) ? ( F0_xn - F0_xnm1 ) / diff : shape_sin(.5*(x+xnm1));
    );
    
    F0_xnm1   = F0_xn;
    xnm1      = x;
  
    antialias
  );

  function fix_dc(x)
  local()
  global()
  instance(DC_fixed, prev)
  (
    DC_fixed=0.999*DC_fixed + x - prev;
    prev=x;
    DC_fixed
  );
  
  function half_delay(x)
  local(nu)
  global()
  instance(y0, d0)
  (
   // Allpass filter to fix phase issues due to the half sample group delay induced by the continuous convolution
   // nu = (1.0-0.5)/(1.0+0.5);
   nu = 0.333333333333333333333333;
   y0 = nu * x + d0 - nu * y0;
   d0 = x;
   y0
  );
  
  function s_delay(x)
  local(nu)
  global()
  instance(y0, d0)
  (
   // Allpass filter to fix phase issues due to the group delay
   nu = 0.036;
   y0 = nu * x + d0 - nu * y0;
   d0 = x;
   y0
  ); 
  
  function shelf_correct(x)
  local()
  global()
  instance(y0, d0, theta, a1, a2, b0, b1)
  (
    y0 = b0 * x + b1 * d0 - a1 * y0;
    d0 = x;
    y0;
  );
  
  function processNonLinearity(inL, inR)
  (
    ( func == 0 ) ? (
      ( slider3 == 2 ) ? (
        outL = ch0.antialiased_tanh_linear(inL);
        outR = ch1.antialiased_tanh_linear(inR);
      ) : ( slider3 == 1 ) ? (
        outL = ch0.antialiased_tanh_rect(inL);
        outR = ch1.antialiased_tanh_rect(inR);
      ) : (
        outL = tanh(inL);
        outR = tanh(inR);
      );
    ) : ( func == 1 ) ? (
      ( slider3 == 0 ) ? (
        outL = clip(inL);
        outR = clip(inR);
      ) : (
        outL = ch0_c.antialiased_clip_rect(inL);
        outR = ch1_c.antialiased_clip_rect(inR);
      );
    ) : ( func == 2 ) ? (
      ( slider3 == 0 ) ? (
        outL = shape_sin(inL);
        outR = shape_sin(inR);
      ) : (
        outL = ch0_s.antialiased_sin_rect(inL);
        outR = ch1_s.antialiased_sin_rect(inR);
      );
    );
  );

  total_gain = preamp * ceiling;
  spl0 *= total_gain;
  spl1 *= total_gain;  
  
  shelf_active = shelf && (oversampling == 1);
  
  slider3 == 1 ? (
    shelf_active ? (
      spl0 = l_delay.s_delay(spl0);
      spl1 = r_delay.s_delay(spl1);
    ) : (
      spl0 = l_delay.half_delay(spl0);
      spl1 = r_delay.half_delay(spl1);
    );
  );    
  
  ( Oversampling == 1 ) ? (
    processNonLinearity(spl0, spl1);
    spl0 = outL;
    spl1 = outR;
  ) : (
    inL = spl0;
    inR = spl1;
    upsampleL.updateUpHist(overSampling, inL);
    upsampleR.updateUpHist(overSampling, inR);
    
    f = 0;
    loop(overSampling,
      f += 1;
      inL = overSampling*upsampleL.upSample(overSampling);
      inR = overSampling*upsampleR.upSample(overSampling);
  
      processNonLinearity(inL, inR);
      downL.updateDownHist(overSampling, outL);
      downR.updateDownHist(overSampling, outR);
      
      ( f == 1 ) ? (
        spl0 = downL.downSample(overSampling);
        spl1 = downR.downSample(overSampling);
      );
    );
  );
  
  slider4 ? (
    spl0 = dc0.fix_dc(spl0);
    spl1 = dc1.fix_dc(spl1);
  );
   
  (lastmode != slider3) ? (
    block = 6;
  );

  (slider3==1) && shelf_active ? (
    spl0 = l_shelf.shelf_correct(spl0);
    spl1 = r_shelf.shelf_correct(spl1);
  );

  ( block > 0 ) ? (
    block = block - 1;
    spl0 = 0;
    spl1 = 0;
  ) : (
    spl0 *= inv_ceiling;
    spl1 *= inv_ceiling;
  );  
  
  lastmode = slider3;
  current_sample += 1;