Fast sine using a UINT32 phase

[neotec]

Today I've reviewed some of my nasm sources and found a nice little routine which calculates the sine value of a UINT32 phase, so basically is does:

Code: Select all

float sinu32(uint32_t phase)
{
 return sinf((float)phase * (3.1415926536f / 2147483648.f))
}

phase step is: 4294967296.f * frequency / FS;

The above code executes 67108864 sine calculations in ~2.9 seconds on my Athlon 64 X2 @ 2.4GHz, my assembler sine does the same job in ~1.3 seconds:

Code: Select all

[bits 32]

[section .bss]

[section .data]
ALIGN	16
fpi:		dd	0x3fc90fdb  ; pi*0.5
frf3:		dd	0xbe2aaaab  ; -1 / 6
frf5:		dd	0x3c088889  ; 1 / 120
frf7:		dd	0xb9500d01  ; -1 / 5040
frf9:		dd	0x3638ef1d  ; 1 / 362880

[section .text]
global _nfastsin2
_nfastsin2:
push    eax
push    ebx
mov     eax, dword[esp + 4 + 8]
mov     ebx, eax
and     ebx, 0x80000000
or      ebx, 0x3f800000
push    ebx
mov     ebx, 0x00800000
bt      eax, 30
sbb     ebx, 0
shr     eax, 7
xor     ebx, eax
or      ebx, 0x3f800000
push    ebx
fld     dword[esp]
fld1
fsubp
fmul    dword[esp + 4]
fmul    dword[fpi]
fld     st0
fmul    st0, st1
fld     st0
fmul    dword[frf9]
fadd    dword[frf7]
fmul    st0, st1
fadd    dword[frf5]
fmul    st0, st1
fadd    dword[frf3]
fmulp   st1, st0
fld1
faddp
fmulp   st1, st0
add     esp, 8
pop     ebx
pop     eax
ret

For using this function you'll need NASM and some header file, like this:

Code: Select all

#ifdef __cplusplus
extern "C"
{
#endif

#ifdef WIN32
	float _cdecl nfastsin2(unsigned __int32 p);
#else
#       include <stdint.h>
#	define _cdecl __attribute__((__cdecl__))
	float _cdecl _nfastsin2(uint32_t p);
#	define nfastsin2(x) _nfastsin2(x)
#endif
#ifdef __cplusplus
}
#endif

Maybe this is useful for somebody
(This routine is the 1-sine version of my SSE additive oscillator core ... I think I posted the sources to the 64 bands add-osc somewhere here in this forum ... but I'll repost it, if somebody is interested)

[duncanparsons]

Just chiming in to bookmark for later

What accuracy is your sin? It looks like some range scaling then a Taylor series - that be about right?

oh, and what are the frf things? is that NASM specific stuff, or have you missed some labels?

I don't recall those other sources you mention, but if you're willing to share it would be interesting to look at

[neotec]

Jup, I forgot to include the taylor constants (post is edited).

About the accuracy: I ran a simple compare to check the maximum error, result is 0.000004f. Should be accurate enough, I think.

About the additive-sine-oscillator, here are the sources. Warning: they aren't thread safe.

The input parameters to these functions are:
phase: integer ... see first post on 'details'
volumes: array of floats, must be aligned to a 16 byte boundary.

About the volumes array: Layout is as follows:
first 4 floats: Level of harmonics 1 - 4
next 4 floats: Volume modulation factor for harmonics 1 - 4
next 4 floats: Level of harmonics 5 - 8
next 4 floats: Volume modulation factor for harmonics 5 - 8
... and so on

Some notes on the code in my first post. This is an optimized assembler version of this C code:

Code: Select all

float sinu32(UINT32 phase)
{
 const float frf3 = -1.0f / 6.0f;
 const float frf5 = 1.0f / 120.0f;
 const float frf7 = -1.0f / 5040.0f;
 const float frf9 = 1.0f / 362880.0f;
 const float f0pi5 = 1.570796327f;
 UINT32 tmp = 0x3f800000 | (phase >> 7);
 if (phase & 0x40000000)
  tmp ^= 0x007fffff;
 float x = (*((float*)&tmp) - 1.0f) * f0pi5;
 float x2 = x * x;
 float asin = ((((frf9 * x2 + frf7) * x2 + frf5) * x2 + frf3) * x2 + 1.0f) * x;
 return (phase & 0x80000000) ? -asin : asin;
}

I just eliminated all conditional branches

[nollock]

All that bit twiddling can be replaced by...

Code: Select all

const double SCALE =  PI / 2147483648.0;

FILD   [phase]
FMUL   [SCALE]
FABS

Which results in [0 .. 80000000 .. FFFFFFFF] mapping to [0 .. PI .. 0]

same as yours does. It just treats the top half as negative and folds up, yours treats it as positive and folds it back down. But the result is the same.

And id bet on it being a fair bit faster.

[neotec]

This would work if the taylor series wasn't only valid for values between 0 and pi/2

So, you have to map 0 - > 1 to
0 -> 0.25 : 0 -> 0.5pi
0.25 -> 0.5 : 0.5pi -> 0
0.5 -> 0.75 : 0 -> 0.5pi (change sign of output)
0.75 -> 1.0 : 0.5pi -> 0 (change sign of output)

So, your solution won't work here.

[valhallasound]

This would work if the taylor series wasn't only valid for values between 0 and pi/2

So, you have to map 0 - > 1 to
0 -> 0.25 : 0 -> 0.5pi
0.25 -> 0.5 : 0.5pi -> 0
0.5 -> 0.75 : 0 -> 0.5pi (change sign of output)
0.75 -> 1.0 : 0.5pi -> 0 (change sign of output)

So, your solution won't work here.

Or you could subtract pi/2 from the results in the previous post. A single floating point addition. The results will be in the range between -pi/2 and pi/2. This works for Taylor series sine approximation - I've done it plenty of times in the past.

Sean Costello

[nollock]

This would work if the taylor series wasn't only valid for values between 0 and pi/2

So, you have to map 0 - > 1 to
0 -> 0.25 : 0 -> 0.5pi
0.25 -> 0.5 : 0.5pi -> 0
0.5 -> 0.75 : 0 -> 0.5pi (change sign of output)
0.75 -> 1.0 : 0.5pi -> 0 (change sign of output)

Ahh I didnt spot that..

But it's easy enough to fix, you just shift the phase left one bit before loading it into the float.

Code: Select all

const double SCALE =  PI / 2147483648.0;

MOV    EAX,[phase]
SHL    EAX,1
MOV    [tmp],EAX
FILD   [tmp]
FMUL   [SCALE]
FABS 

You still need to flip the sign of the output as in your code, obviously.

[neotec]

Ok, I've checked nollocks first post and subtract (as valhallasound stated) pi/2, this new version is about 20% faster.

Code: Select all

;fpip:	dd	0x30c90fdb ; PI / 2^31
global _nfastsin3
_nfastsin3:
fild	dword[esp + 4]
fmul	dword[fpip]
fabs
fsub	dword[fpi]
fld		st0
fmul	st0, st1
fld		st0
fmul	dword[frf9]
fadd	dword[frf7]
fmul	st0, st1
fadd	dword[frf5]
fmul	st0, st1
fadd	dword[frf3]
fmulp	st1, st0
fld1
faddp
fmulp	st1, st0
ret

Results (0x4000000 sin calculations):

Code: Select all

nfastsin2 -> 2097152.000000 : 1.051203 sec
nfastsin3 -> 2097152.000000 : 0.804615 sec
sinf -> 2097152.000000 : 3.566851 sec

[hibrasil]

any chance you can explain how to use this to an asm-noob? Will nasm work on OSX and Windows?

thanks very much,

oli

[neotec]

NASM works on Win/Linux/MacOSX. I'm using the C++ version of eclipse and makefiles for developing C/C++ projects, so to include nasm support my makefiles look like this (win32/nmake):

Code: Select all

asm=nasm.exe
asmflags=-f win32 -Xvc

{..\}.asm.obj:
  @$(asm) $(asmflags) -o $*.obj $<

Then you can add your .asm sources simply by supplying the object name.

The header file to use these routines would look like:

Code: Select all

#ifdef __cplusplus
extern "C"
{
#endif

#ifdef WIN32
	float _cdecl nfastsin2(uint32_t p);
	float _cdecl nfastsin3(uint32_t p);
#else // gcc
#	define _cdecl __attribute__((__cdecl__))
	float _cdecl _nfastsin2(uint32_t p);
	float _cdecl _nfastsin3(uint32_t p);
#	define nfastsin2(x) _nfastsin2(x)
#	define nfastsin3(x) _nfastsin3(x)
#endif

#ifdef __cplusplus
}
#endif

Typedefs:

Code: Select all

#ifdef WIN32
#	define int8_t __int8
#	define uint8_t unsigned __int8
#	define int16_t __int16
#	define uint16_t unsigned __int16
#	define int32_t __int32
#	define uint32_t unsigned __int32
#	define int64_t __int64
#	define uint64_t unsigned __int64
#else
#	include <stdint.h>
#endif