exp(-2 * PI * X) approx

[Christian Budde]

Finally someone agrees. I think I should lead with a post that links to Christian's file.

What is missing in my files is the fact, that there are only a few comments about the final accuracy of the algorithms. I should add typical cases and express the error in dB. So in case of using the fast sin() functions distortion will be added and I need to write the expected amount of distortion (or the highest sidelobe for example).

Unfortunately this takes some time I don't have right now, but maybe someone of you can test this.

Also most of you probably need to convert the code to C(++)...

Christian

[mistertoast]

I have my own sin() approx. I'd like to take spectral snapshots of the various approximations for it.

[Music Engineer]

Hey Christian,

as this thread is about approximations and you are posting here - i have once seen a rational approxmiation for tanh in your implementation of Antti's Moog-filter on musicdsp. i consider this as a really nice approximation and use it myself occasionaly since then because it does not suffer from the problem of unbounded output values (like all polynomials would do). but how did you come up with that formula?

[Music Engineer]

I'm so tempted to start an approximations web page or blog and get opinions and papers from different people. Should I?

a page with a survey of useful approximations would certainly be handy.

[LOSER]

Hey Christian,

as this thread is about approximations and you are posting here - i have once seen a rational approxmiation for tanh in your implementation of Antti's Moog-filter on musicdsp. i consider this as a really nice approximation and use it myself occasionaly since then because it does not suffer from the problem of unbounded output values (like all polynomials would do). but how did you come up with that formula?

I don't know about Christian's approximation but here is how I usually tackle this stuff:

First reduce tanh(x) to an exp(x) function in order to do taylor series (if this is not possible, you need to either do a taylor series directly of the function you try to approximate (but if it can't be be reduced back to an exp(x) function taylor series generally might not perform well), use something like the Pade approximant or fit a polynomial directly) ... so here we go:
f1(x) = tanh(x) =
sinh(x)/cosh(x) =
(exp(x)-exp(-x))/(exp(x)+exp(-x)) =
1-2/( exp(2*x) +1) = (Taylor Series)
1-2/( (1+(2*x)+(4*x^2)/2+(8*x^3)/6+(16*x^4)/24+...) +1) (1)

(1) Now gives as a good looking approximation ... for the positive side ( http://fooplot.com/1-2/((1+(2*x)+(4*x^2)/2+(8*x^3)/6+(16*x^4)/24)+1) ), but nasty for the negative side ... but we can just mirror it later.
So we proceed:

(x^4+2*x^3+3*x^2+3*x)/(x^4+2*x^3+3*x^2+3*x+3) =
which can be written as
H / (H+3)
with H =
(x^4+2*x^3+3*x^2+3*x) = (Horner scheme)
x*(x^3+2*x^2+3*x+3) =
x*(x*(x^2+2*x+3)+3) =
x*(x*(x*(x+2)+3)+3)

Now lets write the C function while remembering we must mirror the postive site to the negative:
=>
float f1(float x)
{
float h = fabs(x);
h = h*(h*(h*(h+2)+3)+3);
if(x>0)
return h/(h+3);
else
return -h/(h+3);
}

This should be a fair approximation. However depending on what is needed I'd usually tune it some more.

For this we'll look at a different example.
E.g. say you only need the range [-1;1] and want as few taylor terms as possible you could get away with (1) reduced to:

f2(x) = 1-2/( (1+(2*x)+(4*x^2)/2+...) +1) (2)

When you now look at (2) ( http://fooplot.com/1-2/((1+(2*x)+(4*x^2)/2)+1) ) and see that it is lower than tanh (in range [-1;1], that's where I most likely try matching the curve by scaling it like this (note by doing so we will loose the asymptotes at 1 and -1 but rather scale them as well, however since this example will be for approximating the range [-1;1] we won't worry about that fact):
f2(x) = (1-2/( (1+(2*x)+(4*x^2)/2+...) +1))*1.1
You really should follow this with a plotter or otherwise it will be hard to follow ... anyway this will cause the bit around 0 to be too steep, so we adjust the curve shape like this:

f2(x) = (1-2/( (1+(1.8*x)+(4*x^2)/2+...) +1))*1.1

A feel for mathematical functions really helps finding what to adjust to get the best results ... now the function is too low around 1 again, so scale the height again:

f2(x) = (1-2/( (1+(1.8*x)+(4*x^2)/2+...) +1))*1.15

Then adjust the steepness again:

f2(x) = (1-2/( (1+(1.7*x)+(4*x^2)/2+...) +1))*1.15

.... then eventually I find that this looks good:

f2(x) = (1-2/( (1+(1.6*x)+(4*x^2)/2) +1))*1.18 (3)

( http://fooplot.com/(1-2/((1+(1.6*x)+(4*x^2)/2)+1))*1.18 )
Transform (3) into:

(1.18*x*(2*x+1.6))/(2*x^2+1.6*x+2) =
(1.18*x*(2*x+1.6))/(x*(2*x+1.6)+2) =
(1.18*T)/(T+2)
with T = x*(2*x+1.6)

Now do the C function which also needs to be mirrored again:
=>
float f2(float x)
{
float t = fabs(x);
t = t*(2*t+1.6);
if(x>0)
return 1.18*t/(t+2);
else
return -1.18*t/(t+2);
}

This way you can do whatever approximation you need, e.g. this:
f3(x) = (x*(2.0*x+1))/(x*(2.0*x+1)+1)
which is minimalistic, preserves the asymptotes at 1 and -1, has a positive deviation from tanh(x) at around 0.5 which is bell shaped, then intersecting tanh(x) at around 0.9 to near 1 at infinity, this would be perfectly suited for a waveshaper

=>
float f3(float x)
{
float t = fabs(x);
t = (t*(2.0*t+1));
if(x>0)
return t/(t+1);
else
return -t/(t+1);
}

P.S. Here is some data about the functions:
TANH APPROXIMATIONS:
Cycles: tanh from std math and f1,2,3 as inline C function
tanh (91360) vs f1 (11328) vs f2 (9359) vs f3 (9859)

Maximum deviation from (std math lib) tanh(x) in range [-1;1]
f1 (0.011594) vs f2 (0.003034) vs f3 (0.039074)

Maximum deviation from (std math lib) tanh(x) in range [-2;2]
f1 (0.021302) vs f2 (0.037185) vs f3 (0.055074)

Maximum deviation from (std math lib) tanh(x) in range [-8;8]
f1 (0.021302) vs f2 (0.163474) vs f3 (0.055074)

Maximum deviation from (std math lib) tanh(x) in range [-100;100]
f1 (0.021302) vs f2 (0.179883) vs f3 (0.055074)


I know there are better approximations and probably Christian's are way better, I just thought I share my method, so whether this approximated functions are any good is up to anyone's personal discretion. Also the measurements (especially the cycle estimation) might not be correct.

[Christian Budde]

I don't know about Christian's approximation but here is how I usually tackle this stuff:

First reduce tanh(x) to an exp(x) function in order to do taylor series (if this is not possible, you need to either do a taylor series directly of the function you try to approximate (but if it can't be be reduced back to an exp(x) function taylor series generally might not perform well), use something like the Pade approximant or fit a polynomial directly) ... so here we go:

f1(x) = tanh(x) =
sinh(x)/cosh(x) =
(exp(x)-exp(-x))/(exp(x)+exp(-x)) =
1-2/( exp(2*x) +1) = (Taylor Series)
1-2/( (1+(2*x)+(4*x^2)/2+(8*x^3)/6+(16*x^4)/24+...) +1) (1)

Bearing in mind, that computers can work perfectly with everything to do with the number of 2, I would definitely stop one step earlier in this approximation!

Code: Select all

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

Now you can easily replace the e^x with a 2^x function (and some scale factors I don't have right now as I am on a vacation).

This 2^x function can be evaluated again using the natural CPU float<->int conversion and a correcting polynomial.

The result is by far closer than any polynomial with a comparable CPU performance (at least of what I found so far).

The polynomials (or rational polynomials as for the other tanh(x) approximation) coefficients have been found using an optimizer. The optimizer was trimmed to either find the minimal error and a continuous error (i.e. no discontinuities in the function). These two sometimes fall together, but not necessarily.

I'll give you some more examples once I'm back from my vacation on tuesday.

Kind regards,

Christian

[antto]

i exported the values of exp(-2.0 * PI * x) where X is 0.0 ... 0.5 into a text file .. at 769 steps/samples
opened it up in CurveExpert
the best fit it found was the "Exponential fit" with error of 0.00000046
Exponential Fit: y=ae^(bx)
a = 0.99999999
b = -6.2831842

anyway, you are searching for something faster, with probably only multiplication and add/subtract ..

i think i like what LOSER posted about the h*(h*(h+2)+3)+3 thing
so i entered a user-defined formula in CurveExpert that looks kind of like that, but let the program decide the constants

Code: Select all

User-Defined Model: y=x*(x*(x+c)+b)+a
Coefficient Data:
a =	0.91647683
b =	-3.9493458
c =	4.1301845

now that's really poor, Error 0.0261..
i added a little bit more "stages" or whatever it's called, to make it a little more complicated, but still with only addition and multiplication:

Code: Select all

User-Defined Model: y=x*(x*(x*(x*(x*(x+f)+e)+d)+c)+b)+a
Coefficient Data:
a =	0.99965458
b =	-6.252406
c =	19.080017
d =	-35.518509
e =	39.252173
f =	-20.146654

now it has twice the number of coefficients, so it'll be twice as expensive?
but the error is: 0.00008737

that's for 0 to 0.5

[LOSER]

Bearing in mind, that computers can work perfectly with everything to do with the number of 2, I would definitely stop one step earlier in this approximation!

Now you can easily replace the e^x with a 2^x function (and some scale factors I don't have right now as I am on a vacation).

This 2^x function can be evaluated again using the natural CPU float<->int conversion and a correcting polynomial.

The result is by far closer than any polynomial with a comparable CPU performance (at least of what I found so far).

The polynomials (or rational polynomials as for the other tanh(x) approximation) coefficients have been found using an optimizer. The optimizer was trimmed to either find the minimal error and a continuous error (i.e. no discontinuities in the function). These two sometimes fall together, but not necessarily.

I'll give you some more examples once I'm back from my vacation on tuesday.

Kind regards,

Christian

Interesting so you're saying that there is a good and fast approximation of 2^x, for x being a float? I guess the trick is in the mentioned correction polynomial ... Anyway can't hardly wait to see that.

[LOSER]

i exported the values of exp(-2.0 * PI * x)

Right, exp(-2*pi*x) ... Sorry I kind of steered this thread off topic ...

now it has twice the number of coefficients, so it'll be twice as expensive?
but the error is: 0.00008737

that's for 0 to 0.5

Well, it all depends on what you got and want, like mystran already hinted it is mostly always better to taylor your approximations for the problem. So e.g. if a multiply-add instruction is available polynomials factored in horner scheme are pretty convenient since you always got a=a*b terms which are single instructions. In case you don't have that something else might be better etc ...

[mistertoast]

"Interesting so you're saying that there is a good and fast approximation of 2^x, for x being a float?"

Yeah, float representation is not very hard to work with using bit masks.

[antto]

true

[LOSER]

Yeah, float representation is not very hard to work with using bit masks.

I know, I was just wondering how this is done for x being a fraction, since for x being integer you can just add to the exponent, but that will give you a stepped curve. However, I see how you can then add some sort of correction to smooth that stepped curve ... which is really interesting and I never thought about that.

[arakula]

[Edit:] Forget it, has already been dealt with on page 2 of this thread. I'd delete this post if I could.

[nollock]

Yeah, float representation is not very hard to work with using bit masks.

I know, I was just wondering how this is done for x being a fraction, since for x being integer you can just add to the exponent, but that will give you a stepped curve. However, I see how you can then add some sort of correction to smooth that stepped curve ... which is really interesting and I never thought about that.

You just split 'x' into integer and fractional components, IE..

x^(a+b) = x^a * x^b

So for 2^x you do this.

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

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

2^f can be done with a polynomial fitted only over the range 0..1

In fact this is how exp() is done on x86 anyway. It's done in two parts, there's an instruction that does 2^x, where x is in the range -1..1, and there's an instruction that does 2^x where x is an integer.

[Christian Budde]

You just split 'x' into integer and fractional components, IE..

x^(a+b) = x^a * x^b

So for 2^x you do this.

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

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

2^f can be done with a polynomial fitted only over the range 0..1

In fact this is how exp() is done on x86 anyway. It's done in two parts, there's an instruction that does 2^x, where x is in the range -1..1, and there's an instruction that does 2^x where x is an integer.

Right. The approximations for 2^x with adjustable accuracy can be found in my pascal unit (mentioned here before).

I guess the x86 solution is done in a very similar manner, but simply with higher accuracy.

Christian