sin/cos implementation in fixed-point

Hello,

Here's my problem. I converted some code from floating-point to fixed
point. Currently, for sin/cos/tan I'm using the floating point math.h
functions. I have two base formats q.22 and q.28 that I use and play
with.

I was initially using 1024 values sin/cos tables (in q.28 format).
However, I'm finding that this is not giving me good-enough precision.
I'm getting abosolute differences of about 0.004-0.008 between the
floating point sin/cos and the lookup tables. The alogrithm (I'm
using) seems very sensitive and I'd like to achieve an absolute error
of 0.000001 (10^-6) correct operaion. For such a precision, the lookup
table seems out of the question because the code is going to run on a
embedded processor with limited cache (8-16K).

So what I'm thinking is maybe to get such a high precision, it may be
possible to write sin/cosine as fixed-point functions (it will be
slower than lookup but hopefully faster than floating point). I was
looking at the implementation of such functions in glibc but couldn't
exactly figure out what's going on in there (mix of assembly/C nasty
stuff).

One possible way is taylor series approximation. However, this won't
work very well since in 32 bits, 10! (ie. 12 terms) will easily
overflow.

I'm sure many of you have come across this problem before in DSP
applications. Since I'm not a DSP guy, I'm not familiar with
solutions. Any of you have any ideas on what I can do? I was reading
somewhere that continued fractions might work. Need to find out more
about that. Any of you have any other approximation methods for
sin/cos/tan that could be done in fixed-point (one can always have as
many base formats as needed) to give the high accuracy that I'm
seeking?

I want to also avoid 64-bit because I plan to then implement parts of
the algorithm in TIE, feed it to Xtensa to generate custom
instructions. However, if 64-bit gives good precision, then I may just
have to stick with it and incur the overhead of doing 64-bit
operations on a 32-bit processor.

I appreciate your advice/help.

pallav

Pallav,

Using a look-up table for sin and cos should work if you interpolate. If 
linear interpolation isn't accurate enough, use quadratic interpolation. 
You need only one table covering a single quadrant for both functions.

A Taylor series need not overflow. Implementing ax^3 + bx^2 + cx + d is 
indeed prone to overflow. It's inefficient, too. Use x(x(ax+b)+c)+d. 
There are fewer multiplications and the small coefficients offset the 
growing x^n.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������
Pallav wrote:

> ... I converted some code from floating-point to fixed point. 

   ...

> I appreciate your advice/help.
> 
> pallav

Try cordic method to generate sin/cos

with regards
praveen

Pallav wrote:

> Hello,
> 
> Here's my problem. I converted some code from floating-point to fixed
> point. Currently, for sin/cos/tan I'm using the floating point math.h
> functions. I have two base formats q.22 and q.28 that I use and play
> with.
> 
> I was initially using 1024 values sin/cos tables (in q.28 format).
> However, I'm finding that this is not giving me good-enough precision.
> I'm getting abosolute differences of about 0.004-0.008 between the
> floating point sin/cos and the lookup tables. The alogrithm (I'm
> using) seems very sensitive and I'd like to achieve an absolute error
> of 0.000001 (10^-6) correct operaion. For such a precision, the lookup
> table seems out of the question because the code is going to run on a
> embedded processor with limited cache (8-16K).
> 
> So what I'm thinking is maybe to get such a high precision, it may be
> possible to write sin/cosine as fixed-point functions (it will be
> slower than lookup but hopefully faster than floating point). I was
> looking at the implementation of such functions in glibc but couldn't
> exactly figure out what's going on in there (mix of assembly/C nasty
> stuff).
> 
> One possible way is taylor series approximation. However, this won't
> work very well since in 32 bits, 10! (ie. 12 terms) will easily
> overflow.
> 
> I'm sure many of you have come across this problem before in DSP
> applications. Since I'm not a DSP guy, I'm not familiar with
> solutions. Any of you have any ideas on what I can do? I was reading
> somewhere that continued fractions might work. Need to find out more
> about that. Any of you have any other approximation methods for
> sin/cos/tan that could be done in fixed-point (one can always have as
> many base formats as needed) to give the high accuracy that I'm
> seeking?
> 
> I want to also avoid 64-bit because I plan to then implement parts of
> the algorithm in TIE, feed it to Xtensa to generate custom
> instructions. However, if 64-bit gives good precision, then I may just
> have to stick with it and incur the overhead of doing 64-bit
> operations on a 32-bit processor.
> 
> I appreciate your advice/help.
> 
> pallav

Analog Devices used to have an excellent example in their 21xx 
literature for 16-bits, but it could be extended.  Basically they used a 
polynomial in x to get 1/4 cycle of a sine function.  They had a 
hard-coded vector of coefficients which they appear to have calculated 
using a LMS fit (I checked, and got the same numbers).  It's 
sorta-quasi-taylor, but the even-valued coefficients do have some weight 
to them which reduces the error over the domain of the function and 
allows you to truncate the series quicker.  They were getting 15-bit 
accuracy with 6 coefficients and about 25 clock cycles (setup, 6x x^2, 
6x MAC).

You should be able to do the same thing with an appropriately expanded 
set of coefficients.  Since you'll be doing all the quadrant-flopping 
anyway you may want to investigate using sin(x - pi/4) -- it may reduce 
the number of coefficients needed.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Jerry Avins wrote:

> Pallav,
> 
> Using a look-up table for sin and cos should work if you interpolate. If 
> linear interpolation isn't accurate enough, use quadratic interpolation. 
> You need only one table covering a single quadrant for both functions.
> 
> A Taylor series need not overflow. Implementing ax^3 + bx^2 + cx + d is 
> indeed prone to overflow. It's inefficient, too. Use x(x(ax+b)+c)+d. 
> There are fewer multiplications and the small coefficients offset the 
> growing x^n.
> 
> Jerry

If you treat the x as Q31 scaled by semicircles then overflow is no 
problem at all.  _Underflow_ may be an issue...

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Tim Wescott wrote:

 > Jerry Avins wrote:
 >
 >> Pallav,
 >>
 >> Using a look-up table for sin and cos should work if you interpolate.
 >> If linear interpolation isn't accurate enough, use quadratic
 >> interpolation. You need only one table covering a single quadrant for
 >> both functions.
 >>
 >> A Taylor series need not overflow. Implementing ax^3 + bx^2 + cx + d
 >> is indeed prone to overflow. It's inefficient, too. Use
 >> x(x(ax+b)+c)+d. There are fewer multiplications and the small
 >> coefficients offset the growing x^n.
 >>
 >> Jerry
 >
 >
 > If you treat the x as Q31 scaled by semicircles then overflow is no
 > problem at all.  _Underflow_ may be an issue...

I've had to go to double precision for one of the multiplies, but that
was with 16-bit words. I think that there's something inherently wrong
with a routine that needs 32-bit trig. It may be easier to write that
way, but I've never seen a case where thinking the problem through
didn't remove the requirement. It roughly amounts to dividing a circle a
mile in diameter into arcs .0001 inches long. The 16-bit code in
http://users.erols.com/jyavins/typek.htm, I treat the numbers as 7.9,
where the 7 includes the sign bit. (Although that code was written to
linearize a thermocouple, it is surprisingly good as a toy sin/cos
routine.)

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

Hello,

Consider having A) a coarsely-spaced table of equally-spaced (cos,sin)
values covering the entire circle, and B) a finely-spaced table of
equally-spaced (cos,sin) values only covering the coarse table
spacing. Then use the equations for cos(small+large), sin(small+large)
in terms of cos(large), sin(large), cos(small), and sin(small) to
compute cos(small+large) and sin(small+large) with 2 multiplies and 1
add each.

With a 100-point 'coarse' table and a 100-point 'fine' table, you have
the equivalent of a 10000-point table (obviously tables do not need to
be equal size).

Few instructions, few operations, little memory, simple addressing,
small numerical error. Additional simplifying tricks possible with a
little thought.

Dirk A. Bell


axl_fugazi@yahoo.com (Pallav) wrote in message news:<baaddb55.0407142104.58ebbd40@posting.google.com>...
> Hello,
> 
> Here's my problem. I converted some code from floating-point to fixed
> point. Currently, for sin/cos/tan I'm using the floating point math.h
> functions. I have two base formats q.22 and q.28 that I use and play
> with.
> 
> I was initially using 1024 values sin/cos tables (in q.28 format).
> However, I'm finding that this is not giving me good-enough precision.
> I'm getting abosolute differences of about 0.004-0.008 between the
> floating point sin/cos and the lookup tables. The alogrithm (I'm
> using) seems very sensitive and I'd like to achieve an absolute error
> of 0.000001 (10^-6) correct operaion. For such a precision, the lookup
> table seems out of the question because the code is going to run on a
> embedded processor with limited cache (8-16K).
> 
> So what I'm thinking is maybe to get such a high precision, it may be
> possible to write sin/cosine as fixed-point functions (it will be
> slower than lookup but hopefully faster than floating point). I was
> looking at the implementation of such functions in glibc but couldn't
> exactly figure out what's going on in there (mix of assembly/C nasty
> stuff).
> 
> One possible way is taylor series approximation. However, this won't
> work very well since in 32 bits, 10! (ie. 12 terms) will easily
> overflow.
> 
> I'm sure many of you have come across this problem before in DSP
> applications. Since I'm not a DSP guy, I'm not familiar with
> solutions. Any of you have any ideas on what I can do? I was reading
> somewhere that continued fractions might work. Need to find out more
> about that. Any of you have any other approximation methods for
> sin/cos/tan that could be done in fixed-point (one can always have as
> many base formats as needed) to give the high accuracy that I'm
> seeking?
> 
> I want to also avoid 64-bit because I plan to then implement parts of
> the algorithm in TIE, feed it to Xtensa to generate custom
> instructions. However, if 64-bit gives good precision, then I may just
> have to stick with it and incur the overhead of doing 64-bit
> operations on a 32-bit processor.
> 
> I appreciate your advice/help.
> 
> pallav

Hi Dirk

Welcome back....haven't noticed your name here in a long time.
Is the absence because you found new work and were ramping up?

Cheers
Bhaskar

"Dirk Bell" <dirkman@erols.com> wrote in message
news:6721a858.0407151847.f13136c@posting.google.com...
> Hello,
>
> Consider having A) a coarsely-spaced table of equally-spaced (cos,sin)
> values covering the entire circle, and B) a finely-spaced table of
> equally-spaced (cos,sin) values only covering the coarse table
> spacing. Then use the equations for cos(small+large), sin(small+large)
> in terms of cos(large), sin(large), cos(small), and sin(small) to
> compute cos(small+large) and sin(small+large) with 2 multiplies and 1
> add each.
>
> With a 100-point 'coarse' table and a 100-point 'fine' table, you have
> the equivalent of a 10000-point table (obviously tables do not need to
> be equal size).
>
> Few instructions, few operations, little memory, simple addressing,
> small numerical error. Additional simplifying tricks possible with a
> little thought.
>
> Dirk A. Bell
>
>
> axl_fugazi@yahoo.com (Pallav) wrote in message
news:<baaddb55.0407142104.58ebbd40@posting.google.com>...
> > Hello,
> >
> > Here's my problem. I converted some code from floating-point to fixed
> > point. Currently, for sin/cos/tan I'm using the floating point math.h
> > functions. I have two base formats q.22 and q.28 that I use and play
> > with.
> >
> > I was initially using 1024 values sin/cos tables (in q.28 format).
> > However, I'm finding that this is not giving me good-enough precision.
> > I'm getting abosolute differences of about 0.004-0.008 between the
> > floating point sin/cos and the lookup tables. The alogrithm (I'm
> > using) seems very sensitive and I'd like to achieve an absolute error
> > of 0.000001 (10^-6) correct operaion. For such a precision, the lookup
> > table seems out of the question because the code is going to run on a
> > embedded processor with limited cache (8-16K).
> >
> > So what I'm thinking is maybe to get such a high precision, it may be
> > possible to write sin/cosine as fixed-point functions (it will be
> > slower than lookup but hopefully faster than floating point). I was
> > looking at the implementation of such functions in glibc but couldn't
> > exactly figure out what's going on in there (mix of assembly/C nasty
> > stuff).
> >
> > One possible way is taylor series approximation. However, this won't
> > work very well since in 32 bits, 10! (ie. 12 terms) will easily
> > overflow.
> >
> > I'm sure many of you have come across this problem before in DSP
> > applications. Since I'm not a DSP guy, I'm not familiar with
> > solutions. Any of you have any ideas on what I can do? I was reading
> > somewhere that continued fractions might work. Need to find out more
> > about that. Any of you have any other approximation methods for
> > sin/cos/tan that could be done in fixed-point (one can always have as
> > many base formats as needed) to give the high accuracy that I'm
> > seeking?
> >
> > I want to also avoid 64-bit because I plan to then implement parts of
> > the algorithm in TIE, feed it to Xtensa to generate custom
> > instructions. However, if 64-bit gives good precision, then I may just
> > have to stick with it and incur the overhead of doing 64-bit
> > operations on a 32-bit processor.
> >
> > I appreciate your advice/help.
> >
> > pallav

> > ... I converted some code from floating-point to fixed point.  
> > I appreciate your advice/help.


Hi Jerry/Praveen/Tim,

Thank you for your responses. I grabbed some papers on cordic and was
able to implement sin/cos (in 32 bits) with the desired accuracy (less
than 10^-6) I was looking for. The algorithm seems to be working fine
now.

Thanks for your time.
pallav

mention not

with regards
praveen