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;

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

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.

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?

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.

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.

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?

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

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!

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.

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.

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).

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.

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?