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

Code: Select all

```
exp2x = exp(2*x)
tanh(x) = (exp2x - 1) / (exp2x + 1)
```

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

(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:

*(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:

*(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)`

Code: Select all

`exp(x) = 1 + x + x^2/(2!) + x^3/(3!) + ...`

Code: Select all

`a1 + a2 * x + a3 * x^2... / b1 + b2 * x + b3 * x^2...`

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;
```

*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:

*(1-Term, no oversampling)*

*(2-Term, no oversampling)*

*(3-Term, no oversampling)*

*(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:

*(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);
```

*(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:

Increasing the order doesn't help much as well:

*(4-Term)*

*(5-Term)*

*(6-Term)*

*(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;
```

*(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:

*(4-Term Continous Error, no oversampling)*

And using a moderate 2x oversampling the alias can be reduced:

*(4-Term Continous Error, 2x oversampling)*

And using a 4x oversampling the alias can be reduced even further:

*(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