View Single Post
Old 10-06-2019, 06:18 AM   #60
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: NL
Posts: 672
Default

I also gave it another shot.

I hunted around for a while for a numerically stable, reasonably fast, but good approximation of the dilogarithm and came across the cern library ROOT:
https://root.cern.ch/root/htmldoc/sr...h.cxx.html#104

Similar to what you did, I had a look for symmetry and noted that aside from a shift, the F1 also has an odd symmetry (although around a non-negative zero point). With this, I did do an implementation of the tanh linear interp, which seems to work alright (though I didn't really put it through its paces enough yet).

Code:
desc:Saike Tanh Saturation with anti aliasing
tags: saturation distortion anti-aliased
version: 1.10
author: Joep Vanlier
changelog: 
  + Re-added block, needed after all.
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.

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)}>Antialias mode
slider4:0<0,1,1>Fix DC?

@slider
preamp      = 10^(slider1/20);
ceiling     = 10^(-slider2/20);
inv_ceiling = 10^(slider2/20);

@sample 
  function F0(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(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 / (diff1*diff1);
      ( xn * ( F0_xn - F0_xnm1 ) - (F1_xn - F1_xnm1) ) * idiff
    ) : (
      .5 * tanh((xn + 2*xnm1)*.33333333333333333333333333333)
    );
    
    term2     = (abs(diff2) > eps) ? (
      idiff = 1 / (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)
  global()
  instance(antialias, F0_xnm1, xnm1,diff)
  (
    (
      F0_xn     = F0(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 fix_dc(x)
  local()
  global()
  instance(DC_fixed, prev)
  (
    DC_fixed=0.999*DC_fixed + x - prev;
    prev=x;
    DC_fixed
  );

  spl0 *= preamp;
  spl1 *= preamp;
  
  spl0 *= ceiling;
  spl1 *= ceiling;  
  
  ( slider3 == 2 ) ? (
    spl0 = ch0.antialiased_tanh_linear(spl0);
    spl1 = ch1.antialiased_tanh_linear(spl1);
  ) : ( slider3 == 1 ) ? (
    spl0 = ch0.antialiased_tanh_rect(spl0);
    spl1 = ch1.antialiased_tanh_rect(spl1);
  ) : (
    spl0 = tanh(spl0);
    spl1 = tanh(spl1);
  );
  
  slider4 ? (
    spl0 = dc0.fix_dc(spl0);
    spl1 = dc1.fix_dc(spl1);
  );
   
  (lastmode != slider3) ? (
    block = 6;
  );

  ( block > 0 ) ? (
    block = block - 1;
    spl0 = 0;
    spl1 = 0;
  ) : (
    spl0 *= inv_ceiling;
    spl1 *= inv_ceiling;
  );

  lastmode = slider3;
There are two things that are a bit meh. One is that when you switch from the other saturators to the linear one, there seems to be a nasty transient at first. Second is that it is a bit more hungry.
__________________
[Tracker Plugin: Thread|Github|Reapack] | [Routing Plugin: Thread|Reapack] | [Filther: Thread|Github|Reapack] | [More JSFX: Thread|Reapack]
sai'ke is offline   Reply With Quote