Old 01-19-2019, 06:01 PM   #41
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: NL
Posts: 1,453
Default

Whoops, yeah, sorry, I chose epsilon too high.

One annoying thing with the method is that it contains a division by a potentially small number (this happens when two subsequent samples have very little difference). To avoid this, you have to replace the solution at that location with what's in the limit there (which can actually be determined analytically). The good part is that the function is relatively well behaved analytically there. The bad part is that it obviously isn't numerically.

Anyways, you have to choose this cutoff point (eps in the code) where you insert the limit solution. I chose it a bit too high.

The second issue is a precision issue and cannot be resolved afaik. But when you're pushing it that hard, it's hardly being used as a normal saturator anymore.

One solution that may be worth investigating is using some sort of lookup table rather than computing the function. Then it can just once be calculated reliably, and numerical issues might not be a thing. Would probably cost some performance tho.

Code:
DO NOT USE THIS VERSION, BETTER VERSION IN LATER POST.

desc:Tanh Saturation with anti aliasing
tags: saturation distortion anti-aliased
version: 1.01
author: Joep Vanlier
changelog: 
  + Tweaked epsilon
  + Tweaked epsilon
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
I have only implemented the rect version, since the linear one depends on Li2 and LUTs aren't so fast in JSFX.

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

@init
bpos=0;

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

@block
blah+=samplesblock;

@sample
spl0=spl0;
spl1=spl1;

@sample 
  function F0(x, em2x)
  local()
  global()
  instance()
  (
    x - log(2/(1 + em2x))
  );
  
  function tanh_prec(x, em2x)
  local() 
  global()
  instance()
  (
    (2/(1+em2x))-1
  );
  
  function tanh(x)
  local()
  global()
  instance()
  (
    (2/(1+exp(-2*x)))-1
  );
  
  function antialiased_tanh_rect(x)
  local(eps, em2x, F0_xn, diff, diffy)
  global(slider4)
  instance(antialias, F0_xnm1, xnm1)
  (
    em2x      = exp(-2*x);
    F0_xn     = F0(x, em2x);
    
    diff      = ( x - xnm1 );
    eps       = 0.000000000000000001;
    antialias = (abs(diff) > eps) ? ( F0_xn - F0_xnm1 ) / diff : tanh_prec(.5*(x+xnm1), em2x);
    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 ? (
    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);
  );
  
  spl0 *= inv_ceiling;
  spl1 *= inv_ceiling;
__________________
[Tracker Plugin: Thread|Github|Reapack] | [Routing Plugin: Thread|Reapack] | [More JSFX: Thread|Descriptions|Reapack]

Last edited by sai'ke; 01-20-2019 at 05:23 AM.
sai'ke is offline   Reply With Quote
Old 01-19-2019, 09:21 PM   #42
JamesPeters
Human being with feelings
 
Join Date: Aug 2011
Location: Near a big lake
Posts: 3,943
Default

Version 1.02 of the tanh saturator starts off with high DC.
JamesPeters is offline   Reply With Quote
Old 01-19-2019, 10:26 PM   #43
ErBird
Human being with feelings
 
Join Date: Jan 2017
Location: Los Angeles
Posts: 1,161
Default

Quote:
Originally Posted by JamesPeters View Post
Version 1.02 of the tanh saturator starts off with high DC.
Indeed. I think this baby needs a little more testing and refinement. Or maybe JS is not the right medium if the linear method yields a better output.

I'm planning to spend some time testing tomorrow.
ErBird is offline   Reply With Quote
Old 01-19-2019, 11:54 PM   #44
JamesPeters
Human being with feelings
 
Join Date: Aug 2011
Location: Near a big lake
Posts: 3,943
Default

Ah well, your original version sounds nice as is aliasing or not, for how I'd use it anyway.
JamesPeters is offline   Reply With Quote
Old 01-20-2019, 01:55 AM   #45
bezusheist
Human being with feelings
 
bezusheist's Avatar
 
Join Date: Nov 2010
Location: Mullet
Posts: 829
Default

see if adding "ext_nodenorm = 1;" @int helps...
__________________
I like turtles
bezusheist is offline   Reply With Quote
Old 01-20-2019, 05:20 AM   #46
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: NL
Posts: 1,453
Default

Code:
desc:Tanh Saturation with anti aliasing
tags: saturation distortion anti-aliased
version: 1.03
author: Joep Vanlier
changelog: 
  + Tweaked epsilon
  + Tweaked epsilon
  + Tweaked epsilon removed factor 0.5
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
I have only implemented the rect version, since the linear one depends on Li2 and LUTs aren't so fast in JSFX.

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

@init
bpos=0;

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

@block
blah+=samplesblock;

@sample
spl0=spl0;
spl1=spl1;

@sample 
  function F0(x, em2x)
  local()
  global()
  instance()
  (
    x - log(2/(1 + em2x))
  );
  
  function tanh_prec(x, em2x)
  local() 
  global()
  instance()
  (
    (2/(1+em2x))-1
  );
  
  function tanh(x)
  local()
  global()
  instance()
  (
    (2/(1+exp(-2*x)))-1
  );
  
  function antialiased_tanh_rect(x)
  local(eps, em2x, F0_xn, diff)
  global(slider4)
  instance(antialias, F0_xnm1, xnm1)
  (
    em2x      = exp(-2*x);
    F0_xn     = F0(x, em2x);
    
    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 ? (
    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);
  );
  
  spl0 *= inv_ceiling;
  spl1 *= inv_ceiling;
I had another look at it with fresh eyes this morning, I was pretty tired last night. There's no reason the numerical issue 'fix' should cause such high frequency spiking even if the epsilon is too high. I noticed that I was missing a factor two somewhere in the before last version, so where I was supposed to replace things near the singularity with an approximation of the function, I was actually putting in samples at half the signal strength. Whoops! The epsilon I was using before was actually better. I think this version is best. No DC when silent, none of that annoying grain when not silent. I tried various tones (high and low frequency) as well as actual audio samples and silence. None of them exhibited pathological behavior as far as I could tell. Let me know if you find another failure mode as I am actually interested in using this architecture in more plugs

Also, I didn't need the ext_nodenorm, but thanks for the heads up. It's a good variable to be aware about
__________________
[Tracker Plugin: Thread|Github|Reapack] | [Routing Plugin: Thread|Reapack] | [More JSFX: Thread|Descriptions|Reapack]

Last edited by sai'ke; 01-20-2019 at 06:13 AM.
sai'ke is offline   Reply With Quote
Old 01-20-2019, 05:38 AM   #47
bezusheist
Human being with feelings
 
bezusheist's Avatar
 
Join Date: Nov 2010
Location: Mullet
Posts: 829
Default

technically there is still DC (+nyquist) in your latest version (-315 dBFS+/-) and adding ext_nodenorm = 1 removes it...fwiw
__________________
I like turtles
bezusheist is offline   Reply With Quote
Old 01-20-2019, 06:20 AM   #48
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: NL
Posts: 1,453
Default

Quote:
Originally Posted by bezusheist View Post
technically there is still DC (+nyquist) in your latest version (-315 dBFS+/-) and adding ext_nodenorm = 1 removes it...fwiw
Yeah, that's not introduced by this specific method though. That's just the denorm noise that you're measuring. I'd rather just keep it on to avoid surprises on the CPU side.
__________________
[Tracker Plugin: Thread|Github|Reapack] | [Routing Plugin: Thread|Reapack] | [More JSFX: Thread|Descriptions|Reapack]

Last edited by sai'ke; 01-20-2019 at 06:25 AM.
sai'ke is offline   Reply With Quote
Old 01-20-2019, 07:33 AM   #49
TonE
Human being with feelings
 
Join Date: Feb 2009
Location: Reaper HAS send control via midi !!!
Posts: 4,031
Default

How could we compare those formula visually, is using R best? How to convert jsfx code to R then? Plotting and zooming to certain interesting areas could help I guess, in comparing and understanding details? Or do we need Matlab/Octave for this?
TonE is offline   Reply With Quote
Old 01-20-2019, 07:52 AM   #50
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: NL
Posts: 1,453
Default

I already looked at it in octave.

With singularity:


As you can see, perpendicular to the singularity, there is not much change. The equation used to replace values near this singularity turns it into this:


The most interesting region is the off center diagonal (where values are still high, but close together). If we zoom in here, we can see that an aggressive cutoff shows better suppression of that singularity without hurting the continuity too much numerically. In the bottom graph I kept the colorbar the same for left and right.


This is why I originally chose that relatively high cutoff. At the time, I didn't realize that I had bug somewhere else tho (it's never where you look)

The first version had a bug. The last one sounds fine afaik.
__________________
[Tracker Plugin: Thread|Github|Reapack] | [Routing Plugin: Thread|Reapack] | [More JSFX: Thread|Descriptions|Reapack]

Last edited by sai'ke; 01-20-2019 at 08:09 AM.
sai'ke is offline   Reply With Quote
Old 01-20-2019, 11:30 AM   #51
ErBird
Human being with feelings
 
Join Date: Jan 2017
Location: Los Angeles
Posts: 1,161
Default

The new version seems very well-behaved. Very nice!

I made a small tweak to allow higher gain levels. It doesn't seem to affect the anti-aliasing ability. Tested up to 90 dB on a 100 Hz sine wave.

Instead of:
Code:
slider3 ? (
  spl0 = ch0.antialiased_tanh_rect(spl0);
  spl1 = ch1.antialiased_tanh_rect(spl1);
) : (
  spl0 = tanh(spl0);
  spl1 = tanh(spl1);
);
I have:
Code:
slider3 ? (
  354 < abs(spl0) ? (
    spl0 = sign(spl0);
  ):(
    spl0 = ch0.antialiased_tanh_rect(spl0);    
  );
  354 < abs(spl1) ? (
    spl1 = sign(spl1);
  ):(
    spl1 = ch1.antialiased_tanh_rect(spl1);    
  );
) : (
  spl0 = tanh(spl0);
  spl1 = tanh(spl1);
);

Last edited by ErBird; 01-20-2019 at 07:16 PM.
ErBird is offline   Reply With Quote
Old 01-20-2019, 04:36 PM   #52
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: NL
Posts: 1,453
Default

Cool! Fix makes sense

Mind if I nick it for the repo?
__________________
[Tracker Plugin: Thread|Github|Reapack] | [Routing Plugin: Thread|Reapack] | [More JSFX: Thread|Descriptions|Reapack]
sai'ke is offline   Reply With Quote
Old 01-20-2019, 05:08 PM   #53
ErBird
Human being with feelings
 
Join Date: Jan 2017
Location: Los Angeles
Posts: 1,161
Default

Quote:
Originally Posted by sai'ke View Post
Cool! Fix makes sense

Mind if I nick it for the repo?
Please do. It's yours. Of course the 354 value can possibly be lowered. I put it just below where the plugin runs out of precision.

Thanks for sharing this anti-aliasing concept with us.

Last edited by ErBird; 01-20-2019 at 07:19 PM.
ErBird is offline   Reply With Quote
Old 01-21-2019, 07:08 AM   #54
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

That anti-alias tanh is interesting... Here is another version of it I came up while playing around with it:

Code:
function tanh_aa(x)
  instance(y)
  local(old_x, old_y, dx)
(
  old_x = this.x;
  old_y = y;

  abs(x) > log(2^511) ? (
    x = sign(x);
    this.x = x * log(2^511);
    y = log(2^511) - log(2);

    x;
  ) : (
    this.x = x;
    y = x - log(2 / (exp(x * -2) + 1));

    dx = x - old_x;
    abs(dx) > 0.0000000001 ? (y - old_y) / dx : 2 / (exp(-(x + old_x)) + 1) - 1;
  );
);
It should be a bit faster than sai'ke's original version (although you will likely not notice, unless you use it many times per sample), and it has ErBird's abs(x)>354 inside the function.

BTW, it seems that this tanh doesn't like sudden large changes in gain, I guess because then the difference between the current and the previous sample becomes too large. A smoother could probably fix this though.

Last edited by Tale; 01-22-2019 at 04:01 AM. Reason: Nitpicked/optimized duplicate sign(x)
Tale is online now   Reply With Quote
Old 01-21-2019, 02:12 PM   #55
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: NL
Posts: 1,453
Default

Cool!

Hey Tale, what did you notice that made you think it doesn't like sudden changes in gain? I've been trying it on pulses of gain, and notice no adverse effects?

I'm curious why you write the value to compare with as log(2^511)? I realize that's close to 354, but is there some reason that would be more performant?
__________________
[Tracker Plugin: Thread|Github|Reapack] | [Routing Plugin: Thread|Reapack] | [More JSFX: Thread|Descriptions|Reapack]
sai'ke is offline   Reply With Quote
Old 01-22-2019, 01:51 AM   #56
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

Quote:
Originally Posted by sai'ke View Post
Hey Tale, what did you notice that made you think it doesn't like sudden changes in gain? I've been trying it on pulses of gain, and notice no adverse effects?
Just send a sine wave (e.g. JS: synthesis/tonegenerator) through the JSFX below, and manually change gain at once between e.g. 0 and 60 dB a few times, and you will see spikes in the output. This doesn't happen with anti-aliasing turned off. I guess it's basically the same issue as with Direct Form biquads.

Code:
desc:Tanh saturation
slider1:0<0,60,1>Gain (dB)
slider2:1<0,1,1{Off,On}>Anti-Aliasing

@init

function tanh(x)
(
  2 / (exp(x * -2) + 1) - 1;
);

function tanh_aa(x)
  instance(y)
  local(old_x, old_y, dx)
(
  old_x = this.x;
  old_y = y;

  abs(x) > log(2^511) ? (
    x = sign(x);
    this.x = x * log(2^511);
    y = log(2^511) - log(2);

    x;
  ) : (
    this.x = x;
    y = x - log(2 / (exp(x * -2) + 1));

    dx = x - old_x;
    abs(dx) > 0.0000000001 ? (y - old_y) / dx : 2 / (exp(-(x + old_x)) + 1) - 1;
  );
);

@slider

a = exp(slider1 * log(10)/20);
b = 1/a;

@sample

x = (spl0 + spl1) * 0.5;
y = aa.tanh_aa(x * a);
slider2 < 0.5 ? y = tanh(x * a);
spl0 = spl1 = y * b;
Quote:
Originally Posted by sai'ke View Post
I'm curious why you write the value to compare with as log(2^511)? I realize that's close to 354, but is there some reason that would be more performant?
No... In fact, using 354 might be safer. But log(2^511) at least hints to where that number came from: The smallest non-denormal 64-bit double floating-point number is 2^(-1022).

BTW, I forgot to say that your one-liner tanh is really great for use in JSFX (and probably C/C++ as well). I was using (exp(2*x)-1)/(exp(2*x)+1), but that requires a variable to store exp(2*x), and it's worse at maintaining precision. So thanks!

Last edited by Tale; 01-22-2019 at 04:04 AM. Reason: Nitpicked/optimized duplicate sign(x)
Tale is online now   Reply With Quote
Old 10-01-2019, 01:24 PM   #57
ErBird
Human being with feelings
 
Join Date: Jan 2017
Location: Los Angeles
Posts: 1,161
Default Update

That 354 comparison and clamping never sat well with me so I recently sought to eliminate it. Above about 60dB the output quickly degraded into aliasing just as bad as a straight tanh.

Stripping down F0, it's clear the loss of precision occurs in e^(-2x) for extreme -x values. For x -> +infinity there is no problem since e^(-2x) -> 0.



F0 is an even function so F0(x) can be changed to F0(abs(x)). This allows you to eliminate the comparison and maintains anti-aliasing up to infinite gain.

Here's the demonstration. Regular tanh is first, then the previous code, then the new one.



The only remaining problem is dealing with the 1/2 sample delay and inherent low-pass effect.

Code:
desc:Tanh AA

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,1,1>Antialias?

@slider

  gain        = 10^(slider1/20);
  ceiling     = 10^(-slider2/20);
  inv_ceiling = 10^(slider2/20);

@sample

  function F0(x)
  (
    abs(x) - log(2/(1 + exp(-2*abs(x))));
  );

  function tanh(x)
  (
    2/(1 + exp(-2*x)) - 1;
  );
  
  function antialiased_tanh_rect(x)
  local(eps, F0_xn)
  instance(xnm1, F0_xnm1, out)
  (
    F0_xn   = F0(x);
    eps     = 0.0000000001;
    out     = eps < abs(x - xnm1) ? (F0_xn - F0_xnm1)/(x - xnm1) : tanh(0.5*(x+xnm1));
    F0_xnm1 = F0_xn;
    xnm1    = x;
    out;
  );

  spl0 *= gain;
  spl1 *= gain;
  
  spl0 *= ceiling;
  spl1 *= ceiling;
  
  slider3 ? (
    spl0 = ch0.antialiased_tanh_rect(spl0);
    spl1 = ch1.antialiased_tanh_rect(spl1);
  ) : (
    spl0 = tanh(spl0);
    spl1 = tanh(spl1);
  );
  
  spl0 *= inv_ceiling;
  spl1 *= inv_ceiling;

Last edited by ErBird; 10-01-2019 at 01:32 PM.
ErBird is offline   Reply With Quote
Old 10-02-2019, 12:29 PM   #58
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: NL
Posts: 1,453
Default

Yep, sounds like a solid solution to the precision problem! Great job

Thanks for sharing it by the way!
__________________
[Tracker Plugin: Thread|Github|Reapack] | [Routing Plugin: Thread|Reapack] | [More JSFX: Thread|Descriptions|Reapack]

Last edited by sai'ke; 10-02-2019 at 12:57 PM. Reason: Edit: Thanks :)
sai'ke is offline   Reply With Quote
Old 10-05-2019, 07:47 PM   #59
ErBird
Human being with feelings
 
Join Date: Jan 2017
Location: Los Angeles
Posts: 1,161
Default

Quote:
Originally Posted by sai'ke View Post
Thanks for sharing it by the way!
No problem. You started the AA ball rolling after all. I wanted to make sure you had the fix.

Anyway, I tried for ages to no avail to get the linear kernel working for tanh. I didn't resort to a lookup table, but tried, with no concern for CPU, to use the power series for Li2 with some abs() and sign() applied accordingly. Alas, it doesn't work.

I did, however, get it working for x/sqrt(1+x^2). Not as nice for audio as tanh, but a good proof that the linear kernel can work. In the end, I think the improvement is not enough to warrant the extra complexity over the rect version. Also it seems more finicky and needs a much larger epsilon and sometimes clicks and pops when fed an already saturated signal. Because of that, I think the latest tanh with AA is the best possible implementation right now.

Code:
desc:Algebraic Saturation with Anti-Aliasing

Author: Erich M. Burg
Adapted from "Tanh Saturation with anti aliasing" JS by Joep Vanlier (Sai'ke)

Uses technique from: "REDUCING THE ALIASING OF NONLINEAR WAVESHAPING USING CONTINUOUS-TIME CONVOLUTION", Parker et al,
Proceedings of the 19th International Conference on Digital Audio Effects (DAFx-16), Brno, Czech Republic, September 5–9, 2016

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 AA,Rectangular Kernel,Linear Kernel}>Mode

@slider

  gain        = 10^(slider1/20);
  ceiling     = 10^(-slider2/20);
  inv_ceiling = 10^(slider2/20);

@sample

  function f(x)
  (
    x/sqrt(1 + x^2);
  );

  function F0(x)
  (
    sqrt(1 + x^2);
  );
  
  function F1(x)
  (
    0.5*(x*sqrt(1 + x^2) - log(sqrt(1 + x^2) + x));
  );

  function antialiased_algebraic_sat_rect(x)
  local(eps, F0_xn, out)
  instance(xnm1, F0_xnm1)
  (
    F0_xn   = F0(x);
    eps     = 0.0000000001;
    out     = eps < abs(x - xnm1) ? (F0_xn - F0_xnm1)/(x - xnm1) : f(0.5*(x+xnm1));
    F0_xnm1 = F0_xn;
    xnm1    = x;
    out;
  );
  
  function antialiased_algebraic_sat_linear(x)
  local(eps, F0_xn, F1_xn, out1, out2, out)
  instance(xnm1, xnm2, F0_xnm1, F0_xnm2, F1_xnm1, F1_xnm2)
  (
    F0_xn   = F0(x);
    F1_xn   = sign(x)*F1(abs(x));
    eps     = 0.0001;
    
    out1    = eps < abs(x - xnm1) ? 
              (
                (x*(F0_xn - F0_xnm1) - (F1_xn - F1_xnm1))/(x - xnm1)^2;
              ):(
                0.5*f((x + 2*xnm1)/3);
              );
    
    out2    = eps < abs(xnm1 - xnm2) ?
              (
                (xnm2*(F0_xnm2 - F0_xnm1) - (F1_xnm2 - F1_xnm1))/(xnm2 - xnm1)^2;
              ):(
                0.5*f((xnm2 + 2*xnm1)/3);
              );
    
    out     = (out1 + out2);
    F0_xnm2 = F0_xnm1;
    F0_xnm1 = F0_xn;
    F1_xnm2 = F1_xnm1;
    F1_xnm1 = F1_xn;
    xnm2    = xnm1;
    xnm1    = x;
    out;
  );

  spl0 *= gain;
  spl1 *= gain;
  
  spl0 *= ceiling;
  spl1 *= ceiling;

  slider3 == 0 ?
  (
    spl0 = spl0/sqrt(1+spl0^2);
    spl1 = spl1/sqrt(1+spl1^2);
  ):
  slider3 == 1 ?
  (
    spl0 = ch0.antialiased_algebraic_sat_rect(spl0);
    spl1 = ch1.antialiased_algebraic_sat_rect(spl1);
  ):
  slider3 == 2 ?
  (
    spl0 = ch0.antialiased_algebraic_sat_linear(spl0);
    spl1 = ch1.antialiased_algebraic_sat_linear(spl1);  
  );
  
  spl0 *= inv_ceiling;
  spl1 *= inv_ceiling;

Last edited by ErBird; 10-05-2019 at 07:53 PM.
ErBird is offline   Reply With Quote
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: 1,453
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] | [More JSFX: Thread|Descriptions|Reapack]
sai'ke is offline   Reply With Quote
Old 10-08-2019, 01:51 PM   #61
solarfall
Human being with feelings
 
Join Date: Sep 2013
Posts: 74
Default

Whoa guys, i didn't receive notifications from this post and thought it was dead long time ago. And now i see 60 replies
Btw, you are far more advanced than me and in the beginning i wasn't even aware of the aliasing problem.
I know some js plugins implement oversampling in their code so it's a metter of copy and paste i guess. Isn't it?
solarfall is offline   Reply With Quote
Old 10-15-2019, 01:41 PM   #62
pepe44
Human being with feelings
 
pepe44's Avatar
 
Join Date: Jul 2013
Location: Portugal
Posts: 1,827
Default

Great thread! always like the saike stuff! Good improvements. Testing your script.
pepe44 is offline   Reply With Quote
Old 10-17-2019, 01:42 AM   #63
solarfall
Human being with feelings
 
Join Date: Sep 2013
Posts: 74
Default

Hello guys,
i did some major corrections and here's a new version. I've reworked the filter which basically is a loPass different for every channel and from left to right.

This plugins works best at 0 VU which corresponds to -18dB. It's supposed to add subtle harmonic distortion so don't feed it with high signal or it will distort badly. I think aliasing should be a minor problem since distiortion is so subtle.

Let me know what you think.

Enjoy!

Code:
desc:AnalogSum

slider1:0<0,2,1{1,2,3}>Saturation Type

@init

a0L = 0.4 + ((rand()/5)-0.05);
a1L = 0.9 + ((rand()/5)-0.05);
a0R = 0.4 + ((rand()/5)-0.05);
a1R = 0.9 + ((rand()/5)-0.05);
dist = 1.6;




function LoPass()
(
  xnL = spl0;
  
  ynL = (a0L * xnL) + (a1L * xn_1L);
  
  xn_1L = xnL;
    
  spl0=ynL;
  
  xnR = spl1;
  
  ynR =  (a0R * xnR) + (a1R * xn_1R);
  
  xn_1R = xnR;
  
  spl1=ynR;
);

function analog_stereo()
(
  m = 0.5*(spl0+spl1);
  s = 0.6*(spl1-spl0);
  spl1 = m + s;
  spl0 = m - s;
);

function add_noise()
(
  spl0 += 0.000005*((rand(1)*2)-1);
  spl1 += 0.000005*((rand(1)*2)-1);
  spl0 += rand(1)/100000*sin(2*$pi*51*t);
  spl1 += rand(1)/100000*sin(2*$pi*51*t);
  spl0 += rand(1)/100000*sin(2*$pi*150*t);
  spl1 += rand(1)/100000*sin(2*$pi*150*t);
  t+=1/srate;
);

function curve1()
(
  spl0 = ((($e^(dist*spl0))-($e^(-dist*spl0)))/(($e^(dist*spl0))+($e^(-dist*spl0))))-(spl0/sqrt(1+(dist*(spl0^2))))+(spl0/dist);
  spl1 = ((($e^(dist*spl1))-($e^(-dist*spl1)))/(($e^(dist*spl1))+($e^(-dist*spl1))))-(spl1/sqrt(1+(dist*(spl1^2))))+(spl1/dist);
  spl0 *= 0.7;
  spl1 *= 0.7;

);

function curve2()
(
  spl0 = sin(tan(spl0))*cos(atan(dist*spl0))*dist*0.9;
  spl1 = sin(tan(spl1))*cos(atan(dist*spl1))*dist*0.9;
  spl0 *= 0.6;
  spl1 *= 0.6;
);

function curve3()
(
  spl0 = (dist*spl0)/(1+abs(2*dist*spl0))*((2+dist)/dist);
  spl1 = (dist*spl1)/(1+abs(2*dist*spl1))*((2+dist)/dist);
  spl0 *= 0.25;
  spl1 *= 0.25;
);



function Saturate ()
(
 
  reg00 == 0 ? curve1();
  reg00 == 1 ? curve2();
  reg00 == 2 ? curve3();

);

@slider
  
  reg00 = slider1; 
  
  
@block

@sample

  add_noise();
  LoPass();
  Saturate();
  analog_stereo();
  
  slider1=reg00;
solarfall is offline   Reply With Quote
Old 10-17-2019, 03:04 AM   #64
pepe44
Human being with feelings
 
pepe44's Avatar
 
Join Date: Jul 2013
Location: Portugal
Posts: 1,827
Default

Nice work. I like the low end on type 1
pepe44 is offline   Reply With Quote
Old 10-19-2019, 03:29 PM   #65
pepe44
Human being with feelings
 
pepe44's Avatar
 
Join Date: Jul 2013
Location: Portugal
Posts: 1,827
Default

What if, we get a version of this but multiband saturator ?
pepe44 is offline   Reply With Quote
Old 10-21-2019, 07:14 AM   #66
Ozman
Human being with feelings
 
Join Date: Feb 2015
Posts: 753
Default

Would aliasing be a factor with a 64bit floating bit depth?
Ozman is offline   Reply With Quote
Old 10-21-2019, 10:49 AM   #67
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: NL
Posts: 1,453
Default

Quote:
Originally Posted by Ozman View Post
Would aliasing be a factor with a 64bit floating bit depth?
Yep. Aliasing happens when the high frequencies no longer 'fit' in the spectrum. The only real solution is a higher samplerate (namely two times the highest frequency in your source).

Draw a fast sinewave, now draw dots on that sinewave at a regular interval that has a lower frequency than your sine wave. You'll see that when the interval is bigger than your fast sinewave, you can connect the dots into a lower frequency sine wave.

So, long story short, not really related to bitdepth, but related to samplerate.
__________________
[Tracker Plugin: Thread|Github|Reapack] | [Routing Plugin: Thread|Reapack] | [More JSFX: Thread|Descriptions|Reapack]
sai'ke is offline   Reply With Quote
Old 10-22-2019, 02:56 PM   #68
pepe44
Human being with feelings
 
pepe44's Avatar
 
Join Date: Jul 2013
Location: Portugal
Posts: 1,827
Default

been testing this JSFX (Analog Summing by solar) and i found an odd stereo shift in the upper frequencies. Anyone can confirm this ?
I did some tweaks on the code but not good results.

Last edited by pepe44; 10-22-2019 at 03:10 PM.
pepe44 is offline   Reply With Quote
Old 10-24-2019, 07:24 AM   #69
Ozman
Human being with feelings
 
Join Date: Feb 2015
Posts: 753
Default

Quote:
Originally Posted by sai'ke View Post
Yep. Aliasing happens when the high frequencies no longer 'fit' in the spectrum. The only real solution is a higher samplerate (namely two times the highest frequency in your source).

Draw a fast sinewave, now draw dots on that sinewave at a regular interval that has a lower frequency than your sine wave. You'll see that when the interval is bigger than your fast sinewave, you can connect the dots into a lower frequency sine wave.

So, long story short, not really related to bitdepth, but related to samplerate.
Can a specified, high samplerate be utilized by the script?
Or is that something that can not be currently done via JSFX?
Ozman is offline   Reply With Quote
Old 10-24-2019, 08:47 AM   #70
solarfall
Human being with feelings
 
Join Date: Sep 2013
Posts: 74
Default

Quote:
Originally Posted by Ozman View Post
Can a specified, high samplerate be utilized by the script?
Or is that something that can not be currently done via JSFX?
Yeah, it can be done by using oversampling. Unfortunately, i'm not as advanced.
solarfall is offline   Reply With Quote
Old 04-09-2021, 06:25 PM   #71
DJ Saint-Hubert
Human being with feelings
 
Join Date: Jan 2013
Posts: 257
Default

Edit: Only the "constant" setting should work, sorry about that. It's because the eps variable is too high for anything to happen with the input volume low as it's supposed to be. Actually it seems like it generates completely different harmonics, i dunno

hope you guys don't mind, I edited the tanh code to use atanh (which is the inverse of tanh, works like an expander instead)

with atanh, the input ranges from -1 to 1 and the output goes from -infinity to +infinity, so lower the initial gain or you will get all sorts of nasty clicks, then the output gain you can raise to somewhere around 30 db to exaggerate the effect (does gain on the difference between the regular and atanh'ed signal)

Code:
desc:Saike Tanh Saturation with anti aliasing, modified by DJ Saint-Hubert
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<-48,12,0.0000001>Gain (dB)
slider2:0<-18,0,1>-Ceiling (dB)
slider3:2<0,2,1{No,Constant,Linear (BETA)}>Antialias mode
slider4:0<0,1,1>Fix DC?
slider5:0<-24,120,0.0000001>Output Gain (dB)

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

@sample 
  function atanh(x)
  local()
  global()
  instance()
  (
    0.5*log((1+x)/(1-x));
  );
  function tanh(x)
  local()
  global()
  instance()
  (
    (0.5*log((1+abs(x))/(1-abs(x))))-abs(x);
  );
  function F0(x)
  local()
  global()
  (
    (x*log(abs(x^2-1)))/(2*abs(x))+x*atanh(abs(x))-(x*abs(x))/2;
  );
  


  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()
  (
    0.5*(abs(x)*log(abs(x^2-1))+(x*log(abs(x+1))-log(abs(x-1))*x)/abs(x)+(-3*log(abs(x)+1)+3*log(1-abs(x))-2*x^2*abs(x))/6+x^2*atanh(abs(x))-abs(x))
  );
  
  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 ) ? (
    spl0a = ch0.antialiased_tanh_linear(spl0);
    spl1a = ch1.antialiased_tanh_linear(spl1);
  ) : ( slider3 == 1 ) ? (
    spl0a = ch0.antialiased_tanh_rect(spl0);
    spl1a = ch1.antialiased_tanh_rect(spl1);
  ) : (
    spl0a = tanh(spl0);
    spl1a = 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 *= ((spl0a*gain)+1);
    spl1 *= ((spl0a*gain)+1);
    spl0 /= preamp;
    spl1 /= preamp;
  );

  lastmode = slider3;
I am doing the 2nd antiderivative right, I hope? It's the antiderivative of the antiderivative I think? I got the formula from a math website, fortunately it didn't involve any dilog functions or anything like that

Would a third antiderivative produce even better results? And how would you interpolate that?

Last edited by DJ Saint-Hubert; 04-10-2021 at 01:14 AM.
DJ Saint-Hubert is offline   Reply With Quote
Old 04-16-2021, 12:10 PM   #72
DJ Saint-Hubert
Human being with feelings
 
Join Date: Jan 2013
Posts: 257
Default

here's the updated script. It only works up to the Constant setting, but with 4x oversampling, the harmonics are almost perfect

https://stash.reaper.fm/v/41798/st-oversampler.jsfx-inc you'll need this too

Code:
desc:Saike Tanh Saturation with anti aliasing, modified by DJ Saint-Hubert
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<-48,12,0.0000001>Gain (dB)
slider2:0<-6,6,0.0000001>-Asym Gain (dB)
slider3:1<0,2,1{No,Constant}>Antialias mode
slider4:0<0,1,1>-Fix DC?
slider5:0<-24,120,0.0000001>Output Gain (dB)
slider6:0<0,1,1{no oversample,4x}>oversample
import st-oversampler.jsfx-inc
@init
pdc_bot_ch = 0; pdc_top_ch = 2;
pdc_delay = 31;
ext_noinit=1;
@slider
preamp      = 10^(slider1/20);
preamp2    = 10^((slider1+slider2)/20);
inv_ceiling = 10^(slider2/20);
gain      = 10^(slider5/20);

@sample 
  function atanh(x)
  local()
  global()
  instance()
  (
    0.5*log((1+x)/(1-x));
  );
  function tanh(x)
  local()
  global()
  instance()
  (
    (0.5*log((1+abs(x))/(1-abs(x))))-abs(x);
  );
  function F0(x)
  local()
  global()
  (
    (x*log(abs(x^2-1)))/(2*abs(x))+x*atanh(abs(x))-(x*abs(x))/2;
  );
  


  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()
  (
    0.5*(abs(x)*log(abs(x^2-1))+(x*log(abs(x+1))-log(abs(x-1))*x)/abs(x)+(-3*log(abs(x)+1)+3*log(1-abs(x))-2*x^2*abs(x))/6+x^2*atanh(abs(x))-abs(x))
  );
  
  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       = .002;
      
    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.0000000000000001;
      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
  );
  
  function process(x) (
                    ( slider3 == 2 ) ? (
                        //spl0a = ch0.antialiased_tanh_linear(spl0);
                        spl0a = ch0.antialiased_tanh_linear(x);
                      ) : ( slider3 == 1 ) ? (
                        //spl0a = ch0.antialiased_tanh_rect(spl0);
                        spl0a = ch0.antialiased_tanh_rect(x);
                      ) : (
                        //spl0a = tanh(spl0);
                        spl0a = tanh(x);
                      );
                      
                      slider4 ? (
                        //spl0 = dc0.fix_dc(spl0);
                        spl0a = dc1.fix_dc(x);
                      );
                       
                      (lastmode != slider3) ? (
                        block = 6;
                      );
                    
                      ( block > 0 ) ? (
                        block = block - 1;
                        //spl0 = 0;
                        x = 0;
                      ) : (
                        //b *= 
                        1;
                      
                       // x /= preamp;
                      );
                     spl0a;
                    );
          
           
           function processr(x) (
                    ( slider3 == 2 ) ? (
                        //spl0a = ch0.antialiased_tanh_linear(spl0);
                        spl1a = ch1.antialiased_tanh_linear(x);
                      ) : ( slider3 == 1 ) ? (
                        //spl0a = ch0.antialiased_tanh_rect(spl0);
                        spl1a = ch1.antialiased_tanh_rect(x);
                      ) : (
                        //spl0a = tanh(spl0);
                        spl1a = tanh(x);
                      );
                      
                      slider4 ? (
                        //spl0 = dc0.fix_dc(spl0);
                        spl1a = dc1.fix_dc(x);
                      );
                       
                      (lastmode != slider3) ? (
                        block = 6;
                      );
                    
                      ( block > 0 ) ? (
                        block = block - 1;
                        //spl0 = 0;
                        x = 0;
                      ) : (
                        //b *= ((spl1a*gain)+1);
                      1;
                        //x /= preamp;
                      );
                      spl1a;
                    );
           
           
           
           

  //spl0 >= 0 ? spl0b = spl0*preamp : spl0b = spl0*preamp2;
  //spl1 >= 0 ? spl1b = spl1*preamp : spl1b = spl1*preamp2;
  spl0b=spl0*preamp;
  spl1b=spl1*preamp;
  slider6 == 0 ? (
  s0=process(spl0b);
  s1=processr(spl1b);
  spl0*=((s0*gain)+1);
  spl1*=((s1*gain)+1);
  ) : (
  s0 = spl0b;
  s1 = spl1b;
  os0.os_up4(s0);
  os1.os_up4(s1);
  
  os0.y3 = process(os0.y3);
  os0.y2 = process(os0.y2);
  os0.y1 = process(os0.y1);
  os0.y0 = process(os0.y0);
  s0 = os0.os_down4();
  
  os1.y3 = processr(os1.y3);
  os1.y2 = processr(os1.y2);
  os1.y1 = processr(os1.y1);
  os1.y0 = processr(os1.y0);
  s1 = os1.os_down4();
  spl0*=((s0*gain)+1);
  spl1*=((s1*gain)+1);
  );
  //spl0 *= ceiling;
  //spl1 *= ceiling;  
  /*
  ( slider3 == 2 ) ? (
    spl0a = ch0.antialiased_tanh_linear(spl0);
    spl1a = ch1.antialiased_tanh_linear(spl1);
  ) : ( slider3 == 1 ) ? (
    spl0a = ch0.antialiased_tanh_rect(spl0);
    spl1a = ch1.antialiased_tanh_rect(spl1);
  ) : (
    spl0a = tanh(spl0);
    spl1a = 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 *= ((spl0a*gain)+1);
    spl1 *= ((spl0a*gain)+1);
    spl0 /= preamp;
    spl1 /= preamp;
  );
*/
  lastmode = slider3;
DJ Saint-Hubert is offline   Reply With Quote
Old 04-17-2021, 03:04 PM   #73
ErBird
Human being with feelings
 
Join Date: Jan 2017
Location: Los Angeles
Posts: 1,161
Default

Hi, DJ. I meant to reply to your first post, sry.

I haven't tried it yet, but looking at the code, I think there are some problems.
What's this function you have as F0? I used Symbolab to find the first and second antiderivatives of arctanh and it's not either of those.
For the constant/rectangular mode you just need the first antiderivative.

Also, you're actually still using tanh when antialias mode is 0:
Code:
spl0a = tanh(spl0);
spl1a = tanh(spl1);
And within antialiased_tanh_rect(x). It should read:
Code:
diff : atanh(.5*(x+xnm1)
I'm curious how this will sound and how it affects dynamics, but aside from that curiosity, I don't think it's going to be a very useful function for a waveshaper.
You usually want to limit dynamic range not dramatically increase it.

There are other types of waveshapers like those in the Geiger Counter guitar pedal which intentionally generate nasty harmonics.
It would be cool to try something like those with this anti-aliasing method.
ErBird is offline   Reply With Quote
Old 04-19-2021, 08:15 AM   #74
DJ Saint-Hubert
Human being with feelings
 
Join Date: Jan 2013
Posts: 257
Default

Quote:
Originally Posted by ErBird View Post
Hi, DJ. I meant to reply to your first post, sry.

I haven't tried it yet, but looking at the code, I think there are some problems.
What's this function you have as F0? I used Symbolab to find the first and second antiderivatives of arctanh and it's not either of those.
For the constant/rectangular mode you just need the first antiderivative.

Also, you're actually still using tanh when antialias mode is 0:
Code:
spl0a = tanh(spl0);
spl1a = tanh(spl1);
And within antialiased_tanh_rect(x). It should read:
Code:
diff : atanh(.5*(x+xnm1)
I'm curious how this will sound and how it affects dynamics, but aside from that curiosity, I don't think it's going to be a very useful function for a waveshaper.
You usually want to limit dynamic range not dramatically increase it.

There are other types of waveshapers like those in the Geiger Counter guitar pedal which intentionally generate nasty harmonics.
It would be cool to try something like those with this anti-aliasing method.
the tanh function is actually atanh, I just thought it was simpler to keep it in. And it's not the antiderivative of atanh(x) because I'm using atanh(abs(x))-abs(x) to get the amount of expansion, which can then be increased if need be

I used this website to find the integral:

https://www.integral-calculator.com/
DJ Saint-Hubert is offline   Reply With Quote
Old 04-19-2021, 01:23 PM   #75
ErBird
Human being with feelings
 
Join Date: Jan 2017
Location: Los Angeles
Posts: 1,161
Default

Oh, interesting. I would rename some functions if tanh is not actually tanh. Also, Li2, F1, antialiased_tanh_linear, etc. are still in the code, just making things more bloated and confusing to read.

What's the reasoning for using atanh(abs(x))-abs(x)? You don't have to answer, I'm just curious what purpose it's serving.

Maybe, alternatively, you want to try sinh(x) if you really want a function that shoots to the stars. Sinh does it a little slower and can accept an input over 0 dB.
ErBird is offline   Reply With Quote
Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off

Forum Jump


All times are GMT -7. The time now is 02:03 AM.


Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2024, vBulletin Solutions Inc.