Tanh approximations

[ToguAudioLine]

It's interesting how you can get close to the shape of the curve but still not match the higher harmonics at all.

Thats really interesting. I would thank Christian Budde for the plots

[Fire Sledge - Ohm Force]

Some times ago I quickly designed a tanh() function approximation :

• tanh(x) ~ sign(x) * (1 - 1 / (1 + |x| + x^2 + a * x^4))

with a = 1 to save a mul (3.0 % max. relative error)
or a = 1.41645 to minimize the maximum relative error (1.7 %)

I haven't checked the harmonics accuracy. Anyway, derivates are continuous up to the third one. It is 12 times faster than the standard tanh() function and almost 30 times faster with SSE packed data. And I'm sure the polynomial can be greatly refined at a small CPU cost (but I'm too lazy to do it).

Code: Select all

double	approx_tanh (double x)
{
	const double	xa = fabs (x);
	const double	x2 = x * x;
	const double	x4 = x2 * x2;
	const double	poly = xa + x2 + x4;
	const double	resa = poly / (1 + poly);
	const double	result = (x < 0) ? -resa : resa;

	return (result);
}

Code: Select all

__m128	approx_tanh (__m128 x)
{
	const __m128	mask_abs = _mm_load_ps (_mask_abs); // 4x 0x7FFFFFFF
	const __m128	one = _mm_load_ps (_one);           // 4x 1.0f

	const __m128	xa = _mm_and_ps (mask_abs, x);
	const __m128	xs = _mm_andnot_ps (mask_abs, x);

	const __m128	x2 = _mm_mul_ps (x, x);
	const __m128	x4 = _mm_mul_ps (x2, x2);

	const __m128	poly = _mm_add_ps (_mm_add_ps (xa, x2), x4);
	const __m128	resa = _mm_div_ps (poly, _mm_add_ps (one, poly));
	const __m128	result = _mm_or_ps (resa, xs);

	return (result);
}

[Christian Budde]

Hi Fire Sledge,

the function generates higher harmonics than necessary for an original tanh. This probably leads to more alias then necessary and such. Comparing both versions the version with 1 instead of 1.41645 performs better from a harmonic point of view.

It's pretty fast as all the approximations based on a rational polynomial function, but not that accurate.

Christian

[Fire Sledge - Ohm Force]

I did some spectrum measurements this morning and came with a better approximation based on the form I proposed yesterday.

The previous function looked like this:



This is the spectrum of a 2*sin(t) signal passed through the function. The green curve is the original tanh(), and the blue one the approximation.

So I changed the function like this :

• tanh(x) ~ sign(x) * (1 - 1 / (1 + |x| + x^2 + 0.66422417311781 * |x|^3 + 0.36483285408241 * x^4))

giving the following results:


The shape in the time domain suffered a bit but the extra-harmonics are better than the first one (same slope, but more than 30 dB of improvement).

EDIT : typo in the formula

[antto]

I did some spectrum measurements this morning and came with a better approximation based on the form I proposed yesterday.

The previous function looked like this:



This is the spectrum of a 2*sin(t) signal passed through the function. The green curve is the original tanh(), and the blue one the approximation.

So I changed the function like this :

• tanh(x) ~ sign(x) * (1 - 1 / (1 + |x| + x^2 + 0.66422417311781 * |x|^3 + 0.36483285408241 * x^4))

giving the following results:


The shape in the time domain suffered a bit but the extra-harmonics are better than the first one (same slope, but more than 30 dB of improvement).

EDIT : typo in the formula

hi Fire Sledge, i entered some tanh() data into CurveExpert
i used only the positive part, from 0 to where it hits 1.0
used your approximitation as a start, and added one more "complication there"
..left CurveExpert to find the perfect coefficients
here's my result:

User-Defined Model: y=(1 - 1 / (1 + x + x^2 + a * x^3 + b * x^4 + c * x^7))
Coefficient Data:
a = 0.58576695
b = 0.55442112
c = 0.057481508

so that's the positive part up to 1.0
the sign(x) and abs(x) things should be used
the input might need clipping so that it doesn't "escape" out of the approximitation.. blah
matches quite close, i can't test the CPU Usage..

[antto]

and a pseudo function..
double tanh_approx2(double x)
{
double xa = abs(x);
double x2 = xa * xa;
double x3 = xa * x2;
double x4 = x2 * x2;
double x7 = x3 * x4;
double res = (1.0 - 1.0 / (1.0 + xa + x2 + 0.58576695 * x3 + 0.55442112 * x4 + 0.057481508 * x7));
return (x > 0 ? res : -res);
}