Tanh approximations

DSP, Plug-in and Host development discussion.
KVRAF
1538 posts since 14 May, 2004 from Europe

Post Wed Sep 09, 2009 9:46 am

Hi,

based on the discussion in this thread, I decided to do some comparison of different tanh() approximations I came up with some time ago for my open source project.

MOTIVATION
A tanh() can often been found when trying to emulate a small signal schematics. Several interesting papers like the moog ladder filter by Antti uses the tanh() in several places and it is often suggested to reduce the amount of calculations by using a polynomial approximation, without loosing a word how this approximation does look like.

SOME MATH
As said in the other thread, a tanh() can be solved using exponential functions. For more information see the wikipedia article about hyperbolic functions

Now here's our starting point:

Code: Select all

tanh(x) = (exp(2*x) - 1) / (exp(2*x) + 1)
The exp(2X) obviously only needs to be calculated once:

Code: Select all

exp2x = exp(2*x)
tanh(x) = (exp2x - 1) / (exp2x + 1)
Which is at least already more efficient then the default Delphi implementation.

LOSSY CODE OPTIMIZATION
Now when we look at this piece of code we can see that there are 2 'ugly' parts. First the exp() function uses quite a huge number of cycles (about 100) and the division (about 10 cycles). Since the exp() function is dominant it's important to tweak our algorithm there.

The exp(x) function can be translated in 2^(c * x) with c = 1 / ln(2). This c can also contain any other factor we need. In the above case we need to calculate exp(2*x), so we can use a factor c' = 2 / ln(2) to simplify our function even further.

Now (as written in the other thread)
nollock wrote:for 2^x you do this:

i = Floor(x)
f = x-i;
result = 2^i * 2^f;

where 2^i can be done with bit masking or bit shifting.

2^f can be done with a polynomial fitted only over the range 0..1
IMPLEMENTATION DETAILS
There was a post by Laurent de Soras some time ago: http://www.musicdsp.org/showone.php?id=63

I wrote a small paper about how to extend his little algorithm towards different levels of accuracy, but I lost that paper due to a head crash of my HDD.

However, the code can still be found in my open source project in this unit http://delphiasiovst.svn.sourceforge.ne ... iew=markup

The problem with the polynomial for 2^f lies in the fact that the value for 0 and 1 does not necessary be the same value and when you glue each section together you get discontinuities.

The polynomial used by Laurent de Soras had the advantage to be not only continuous, but also having the first derivative being continuous. This reduces the amount of extra harmonics being created, but we come to this later.

Using this approximation for the exp() function it is easy to find a series of tanh() approximations. The accuracy directly depends on the order of the polynomial for the 2^f part. The higher the better, but also more CPU usage. At some point you are on par with the x86 implementation (even though you might not necessarily be as close as that).

SOME PSYCHOACOUSTICS
Though in most cases you don't really need the accuracy the x86 implementation provides. If it is inaudible or even falls below the threshold of what can be possible be represented by the DA converters you're probably already too close.

To evaluate what is audible and what is not you can either perform listening tests or estimate by some simple rules. First try to avoid aliasing like plaque. This artifact creates inharmonic frequencies, that only depend on the samplerate (and nothing natural). It does sound unnatural in all cases and there is no excuse. Unfortunately you can't do it completely alias-free, but you can try to make it inaudible. In case of the tanh() function you could try to process a 10 kHz signal looking at the components below that frequency.

Image
(10 kHz excitation [green] & alias frequencies)

This alias can be avoided by oversampling and filtering (aka band limiting).

Now let's have a look at something like output = tanh(2 * input) with input being a 1 kHz sine (NOTE: the factor 2 is used to generate a larger amount of harmonics).

SOME MEASUREMENT RESULTS
The first plot of the pure tanh() function doesn't contain much aliasing mainly because the excitation frequency is quite low compared to the sampling rate:

Image
(Pure digital Tanh(2x) magnitude spectrum without oversamling

Some aliasing can easily be spotted by the intermodulation products around the 1 kHz. As mentioned before the alias will increase for higher gain factors and higher frequencies. It can only be avoided by proper oversampling and filtering.

The oversampling filters in these examples can be adjusted between a factor of 1x-16x using 16th order butterworth filters and a transition bandwidth of 2000 Hz. In case of a typical sampling rate of 44.1 kHz this means an attenuation of -3dB at the cutoff point (about 20 kHz here). Even though the butterworth filters do not produce ripple on their own, a slight ripple around 20 kHz can be detected in this application (with a tanh() waveshaper).

Finally with 16x oversampling the most alias artifacts can be removed:
Image
(16x oversampled tanh(2x) function using a 16th order butterworth oversampling filter)

EVALUATION OF APPROXIMATION
Now let's have a look at some typical rational polynomial approximations.

As explained in detail in the Music-DSP archive the series is motivated by the same basic exponential series:

Code: Select all

tanh(x) = (exp(2x) - 1) /(exp(2x) + 1)
The exp functions are now exchanged by the first values of it's taylor series:

Code: Select all

exp(x) = 1 + x + x^2/(2!) + x^3/(3!) + ...
This results directly in a rational polynomial approximation, which have the general form of

Code: Select all

a1 + a2 * x + a3 * x^2... / b1 + b2 * x + b3 * x^2...
However since the same exp(2x) function is used in the nominator and denominator we can reduce the formula a little bit.

Here are the basic formulas for the first 4 orders/terms (PASCAL code):

Code: Select all

function FastTanh2Like1Term(const Input: Single): Single;
begin
 Result := Input / (abs(Input) + 3);
end;

function FastTanh2Like2Term(const Input: Single): Single;
var
  a, b: Single;
begin
 a := abs(Input);
 b := 3 + a;
 Result := (Input * b) / (a * b + 6);
end;

function FastTanh2Like3Term(const Input: Single): Single;
var
  a, b: Single;
begin
 a := abs(Input);
 b := (6 + a * (3 + a));
 Result := (Input * b) / (a * b + 12);
end;

function FastTanh2Like4Term(const Input: Single): Single;
var
  a, b: Single;
begin
 a := abs(Input);
 b := 12 + a * (6 + a * (3 + a));
 Result := (Input * b) / (a * b + 24);
end;
Due to the fact that the series has been truncated an error is made which leads in general to a smaller absolute function values. This can be seen quite good in the time domain:

Image
Time-Domain display processing a 1 kHz sine wave

The red curve is the first order/term version while the blue curve is the fourth order/term version. The suggested 2nd order version in the post at music-dsp does not perform well here as its values are quite far off the original tanh() version (even though the harmonic content seems to fit as can be seen below)

But now let's have a look at the harmonic distortion created by this waveshaper:

Image
(1-Term, no oversampling)

Image
(2-Term, no oversampling)

Image
(3-Term, no oversampling)

Image
(4-Term, no oversampling)

From the frequency response it seems as if we get an improvement only for every odd order (1, 3, 5...), while there is an improvement in the time domain as could be seen above.

Now let's have a look at the 3rd term/order approximation:
Image
(3-Term, 2x oversampling)

Bearing in mind that this approximation is not too close to a real tanh() in the time domain, it already gives a quite good approximation of the harmonic content.

FITTED RATIONAL POLYNOMIAL APPROXIMATION
Based on the 3rd order/term code, the coefficients can be changed to match the tanh() function in the time domain. This can be necessary as the approximation function has a systematic error of reproducing too small absolute values.

The coefficients below has been found using an optimizer to fit the tanh function as close as possible to the rational polynomial expression over a range of x = -5..+5. A more narrow range could be evaluated, but the results are similar.

Now let's consider the following simple 3-term equation:

Code: Select all

a = abs(Input);
b = 1.26175667589988239 + a *
    (-0.54699348440059470 + a *
    ( 2.66559097474027817));
Result = (b * Input) / (b * a + 1);
Image
(3-Term approximation using optimized coefficients)

We can see that this generates a lot of harmonics not present in the original tanh() function nor the presented approximations above. If not oversampled, these harmonics do lead to aliasing.

But even if oversampled 16x(!) the alias can't be vanished completely:
Image

Increasing the order doesn't help much as well:

Image
(4-Term)

Image
(5-Term)

Image
(6-Term)

Image
(7-Term)

In contrast to this we can look at the proposed approximation using the floating point format. It uses a function presented in the other thread as well:

Code: Select all

const 
  CP2ContError3 : array [0..2] of Single = ( 
    6.93282526441610814E-1, 
    2.42201488582370950E-1, 
    5.50043626970249666E-2); 

function FastPower2ContinousError3(Value: Single): Single; 
var 
  IntCast : Integer absolute result; 
begin 
 IntCast := round(Value); 
 Value := Value - IntCast; 
 IntCast := ($7F + Intcast) shl 23; 
 Result := Result * (1 + 
            Value * (CP2ContError3[0] + 
            Value * (CP2ContError3[1] + 
            Value *  CP2ContError3[2]))); 
end;
The plot looks like this
Image
(3-Term Continous Error)

We can still see some aliasing, but already lower than with the rational polynomial.

This is even better with a 4-Term solution:

Image
(4-Term Continous Error, no oversampling)

And using a moderate 2x oversampling the alias can be reduced:
Image
(4-Term Continous Error, 2x oversampling)

And using a 4x oversampling the alias can be reduced even further:
Image
(4-Term Continous Error, 4x oversampling)

It seems as if there is still some kind of unsymmetrical behaviour as the first harmonic is excited as well, but if we ignore this fact it does already look quite good and still way faster than the original tanh() function provided by most math libraries.

TIMINGS
If we look at the timings I can only provide some simple measurements here. I processed several seconds of random noise with the waveshaper and based on that I calculate the lowest ratio of CPU time vs. audio time [in percent]

So the following values should roughly give an idea how fast the proposed functions are:

Delphi Math Tanh() = 6.9% CPU usage
Mathematical Optimized Tanh() = 2.0% CPU usage
TanhOpt3..7() = 1.0% CPU usage (mainly ruled by division)
TanhContinousError2() = 0.8% CPU usage
TanhContinousError7() = 1.1% CPU usage
TanhLike1() = 0.4% CPU usage
TanhLike4() = 0.7% CPU usage

I still have to provide more exact timings, but not today anymore.

I hope you can make use of this post,

Christian
Last edited by Christian Budde on Fri Sep 11, 2009 1:15 am, edited 8 times in total.

KVRian
1339 posts since 26 Aug, 2005 from Netherlands

Post Wed Sep 09, 2009 10:09 am

Still reading.. but here's already a huge :tu:

KVRAF
3404 posts since 15 Sep, 2002

Post Wed Sep 09, 2009 1:20 pm

Really, really cool.
Swing is the difference between a drum machine and a sex machine.

KVRAF

Topic Starter

1538 posts since 14 May, 2004 from Europe

Post Thu Sep 10, 2009 12:16 pm

Today I replaced some graphics, but there are still some explanations missing. I'll add them as soon as I am able to.

Christian

KVRAF
8384 posts since 11 Apr, 2003 from now on the flat

Post Thu Sep 10, 2009 3:23 pm

Christian - tour de force.. destined to become a classic posting, just fantastic. Good stuff :)
Image

KVRAF

Topic Starter

1538 posts since 14 May, 2004 from Europe

Post Fri Sep 11, 2009 1:21 am

Today I still added some more stuff, but there are still some information missing. Maybe tomorrow I'll have the post entirely done...

KVRAF
3404 posts since 15 Sep, 2002

Post Sat Sep 12, 2009 9:19 am

I'm trying to figure out your methods. Could you walk through how you'd evaluate this approximation for tanh?

Y=X/(A+B*X*X)+C*X

A = 0.11605589854E+01
B = 0.37011709966E+00
C = 0.10710006710E+00

What tools and programs are you using? Given a function, how are you making your frequency plots?

In other words, can you teach me to fish?
Swing is the difference between a drum machine and a sex machine.

KVRAF

Topic Starter

1538 posts since 14 May, 2004 from Europe

Post Sun Sep 13, 2009 3:14 am

mistertoast wrote:I'm trying to figure out your methods. Could you walk through how you'd evaluate this approximation for tanh?
Most approximations have been evaluated by minimizing the maximum error between a tanh() function and its approximation. However, care must be taken that you either use the maximum absolute error or the maximum percent error.

The optimizer used for this isn't really important. It should only be able to find the global minimum rather than a local minimum.

For my evaluations I've been using an evolution algorithm called 'Differential Evolution'. You should find a whole lot of information about that by googling that name.
mistertoast wrote:What tools and programs are you using? Given a function, how are you making your frequency plots?
The tool for plotting is my very own VST Plugin Analyser running a THD test and using the SVG export. For everything else I am solely using the code of my open source project 'Delphi ASIO & VST Project'.

The main thing though that makes these evaluations easy is simply experience. I've tested so many approximations in the past years, that I am quite used to what tricks to use and what to look for.

Kind regards,

Christian

KVRian
1153 posts since 10 Dec, 2003

Post Sun Sep 13, 2009 4:51 pm

Christian Budde wrote: Most approximations have been evaluated by minimizing the maximum error between a tanh() function and its approximation. However, care must be taken that you either use the maximum absolute error or the maximum percent error.
If you're doing piecewise polynomials then you may also want to add extra cost onto the end points and the gradients of the end points. For example...

y = 2^x; // over range 0..1

diferentiate to get this..

dy/dx = ln(2)*2^x;

so when

x = 0, y = 1, dy/dx = ln(2)*1;
x = 1, y = 2, dy/dx = ln(2)*2;

If you add a cost multiplier of 1000 or even 100,000 times for x=0, and x=1, and for the grandient at those points and you'll almost eliminate distortion from the discontinuity between polynomial segments. Well you can basically tweak the cost multiplier until the discontinuity has less effect than other distortion.

KVRian
1296 posts since 12 Apr, 2002

Post Mon Sep 14, 2009 5:14 am

nollock wrote:If you add a cost multiplier of 1000 or even 100,000 times for x=0, and x=1, and for the grandient at those points and you'll almost eliminate distortion from the discontinuity between polynomial segments. Well you can basically tweak the cost multiplier until the discontinuity has less effect than other distortion.
Except that there are still discontinuities of higher orders...

Regards,
{Z}

KVRian
1153 posts since 10 Dec, 2003

Post Mon Sep 14, 2009 8:17 am

Z1202 wrote:
nollock wrote:If you add a cost multiplier of 1000 or even 100,000 times for x=0, and x=1, and for the grandient at those points and you'll almost eliminate distortion from the discontinuity between polynomial segments. Well you can basically tweak the cost multiplier until the discontinuity has less effect than other distortion.
Except that there are still discontinuities of higher orders...
In my experience that doesnt matter. Those higher order discontinuities are likely to have no more effect than the distortion from a less than perfect fit anyway.

KVRian
1153 posts since 10 Dec, 2003

Post Tue Sep 15, 2009 9:44 pm

Ok FWIW i dug out the code and fitted a polynomial to 2^x over the range 0..1, this has 1000000 the cost on the end points, and the gradients of the end points, so it should kill distortion from the inter segment discontinuities.

0.079441542*x^3 + 0.22741128*x^2 + 0.69314718*x + 1.0000000

KVRAF

Topic Starter

1538 posts since 14 May, 2004 from Europe

Post Tue Sep 15, 2009 11:14 pm

Hi there,

out of curiosity I performed the same tests on your set of coefficients and the result was far worse:

Image

I am sure that I have done it right so far, as the shape itself for the resulting function is quite close to a tanh(). There might be some effects due to coefficient truncation, especially for the x^3 coefficient, but I doubt that it will change much.

Even with oversampling it won't get better. Maybe I've done something entirely wrong, but I doubt it.

Perhaps you could post the complete function (without coefficient truncation), just to be sure.

Christian

KVRAF
3404 posts since 15 Sep, 2002

Post Wed Sep 16, 2009 5:46 am

Can you chart my function, Christian?
Swing is the difference between a drum machine and a sex machine.

KVRAF

Topic Starter

1538 posts since 14 May, 2004 from Europe

Post Wed Sep 16, 2009 6:07 am

Hi Mr. Toast,

sorry I missed your function. It looks quite good. It is quite close to the original tanh-function, while having a quite similar harmonic structure.

Find the plot here:
Image

However, the harmonics do not match exactly. They decay more than the original tanh() function, which would explain the reduced aliasing. In the above example the 8 kHz harmonic is already 1 dB too low (which can't be seen in the plot obviously). Near the samplerate it's far more off.

However, it's a nice function.

Christian

Return to “DSP and Plug-in Development”