Reply by Vladimir Vassilevsky●February 25, 20102010-02-25
dewagter wrote:
> Here is a possible fixed point version of the above code. The result is
> accurate within +/-1 degree.
:))))))) +/- 1 km
> If +/-5 degree errors are acceptable then the
> r^3 term can be left out with coeff_1b becoming -45. This function returns
> degrees and not radians.
The best way to approximate atan2() is LUT with trivial linear
interpolation. It takes LUT as small as 32 elements to get the 16-bit
accurate result.
Vladimir Vassilevsky
DSP and Mixed Signal Design Consultant
http://www.abvolt.com
Reply by dewagter●February 25, 20102010-02-25
Here is a possible fixed point version of the above code. The result is
accurate within +/-1 degree. If +/-5 degree errors are acceptable then the
r^3 term can be left out with coeff_1b becoming -45. This function returns
degrees and not radians.
-------------------------------------------------------------------
#define MULTIPLY_FP_RESOLUTION_BITS 15
int16_t atan2_fp(int16_t y_fp, int16_t x_fp)
{
int32_t coeff_1 = 45;
int32_t coeff_1b = -56; // 56.24;
int32_t coeff_1c = 11; // 11.25
int16_t coeff_2 = 135;
int16_t angle = 0;
int32_t r;
int32_t r3;
int16_t y_abs_fp = y_fp;
if (y_abs_fp < 0)
y_abs_fp = -y_abs_fp;
if (y_fp == 0)
{
if (x_fp >= 0)
{
angle = 0;
}
else
{
angle = 180;
}
}
else if (x_fp >= 0)
{
r = (((int32_t)(x_fp - y_abs_fp)) << MULTIPLY_FP_RESOLUTION_BITS) /
((int32_t)(x_fp + y_abs_fp));
r3 = r * r;
r3 = r3 >> MULTIPLY_FP_RESOLUTION_BITS;
r3 *= r;
r3 = r3 >> MULTIPLY_FP_RESOLUTION_BITS;
r3 *= coeff_1c;
angle = (int16_t) ( coeff_1 + ((coeff_1b * r + r3) >>
MULTIPLY_FP_RESOLUTION_BITS) );
}
else
{
r = (((int32_t)(x_fp + y_abs_fp)) << MULTIPLY_FP_RESOLUTION_BITS) /
((int32_t)(y_abs_fp - x_fp));
r3 = r * r;
r3 = r3 >> MULTIPLY_FP_RESOLUTION_BITS;
r3 *= r;
r3 = r3 >> MULTIPLY_FP_RESOLUTION_BITS;
r3 *= coeff_1c;
angle = coeff_2 + ((int16_t) (((coeff_1b * r + r3) >>
MULTIPLY_FP_RESOLUTION_BITS)) );
}
if (y_fp < 0)
return (-angle); // negate if in quad III or IV
else
return (angle);
}
Reply by robert bristow-johnson●April 15, 20052005-04-15
in article k6ednRwgMo1ucsLfRVn-3w@rcn.net, Jerry Avins at jya@ieee.org wrote
on 04/15/2005 12:46:
> Where in your code to you get the numeric value of atan()?
just thought i might offer a little tidbit here for your toolbox:
i picked up a little on Rick Lyons' trick on pp. 547-549 and took part of it
a few further steps. in it he says "Unfortunately, because it is such a
non-linear function, the arctangent is resistant to accurate reasonable-
length polynomial approximations." and while that is true, the denominator
function, g(x), of this:
arctan(y) = y/g(y) -1 <= y <= 1
lends itself quite nicely to the remez exchange algorithm to get as good of
an approximation as you're willing to pay for (given that |y| <= 1). Rick
had a "first" order approximation:
arctan(y) ~= y/(1 + 0.28125*y^2) -1 <= y <= 1
i just wanted to take the denominator a little further and get a better
approximation. since arctan(x) is an odd symmetry function, g(x) has to be
even symmetry and thus has only even powers of x. so the polynomial to fit
is:
f(x) = g(sqrt(x)) = sqrt(x)/arctan(sqrt(x)) 0 < x <= 1
and a really good 4th order approximation i got was:
f(x) ~= 1.0
+ 0.33288950512027 * x
+ -0.08467922817644 * x^2
+ 0.03252232640125 * x^3
+ -0.00749305860992 * x^4
(note that my 0.33289 is a little bigger than Rick's 0.28125 .)
of course you still have to deal with the cost of division (actually, you
have two of them, one of which Rick was able to avoid by limiting his order
to y^2) and twisting the math a little to put it in the correct "octant" to
do atan2().
might be useful. perhaps not.
--
r b-j rbj@audioimagination.com
"Imagination is more important than knowledge."
Reply by pea●April 15, 20052005-04-15
>atan2(u,v) returns angles from 0 to (but not including) 360 degrees,
>whereas atan(u/v) returns angles from -90 to +90 degrees. With either
>function, it is necessary to compute an arctangent. I don't see one in
>your code. One way to implement atan2() is getting atan(abs(u).abs(v)),
>then translating the result to the proper quadrant determined by sgn(u)
>and sgn(v). Where in your code to you get the numeric value of atan()?
>
>Jerry
Hi Jerry,
My code is copied and pasted directly from the example here:
http://www.dspguru.com/comp.dsp/tricks/alg/fxdatan2.htm
Which is by that "original author" above.
The page title says that this computes atan2 (Fast fixed-pt. atan2
calculation with self normalization) and the description/content seems to
verify.
The quadrant is determined and factored in within the equation already
(hence the branches [ if (x>=0) ] and [ if (y<0) ].
Or am I missing something and this is only part of the code?
This message was sent using the Comp.DSP web interface on
www.DSPRelated.com
Reply by Jerry Avins●April 15, 20052005-04-15
pea wrote:
> Hello all,
>
> Sorry to drag up old posts, but I found the original link and have been
> trying to implement atan2 in fixed point. AFter it not working, I found
> this post about the same algorithm and thought it a good place to start.
>
> I'm using 16.16 fixed point. My version of the atan2 from the original
> link follows. Please note - FIX90 and FIX270 are defines for PI/4 and
> 3*PI/4 in degrees. Looking at the code I thought it didn't matter that I
> use degrees rather than radians, as I am expecting the result in degrees
> anyway. Am I correct?
>
> fixed fixatan2(fixed y, fixed x){
> fixed abs_y = y;
> fixed r, angle;
>
> if (x==0){
> if (y==0) return 0L;
> else return ((y < 0) ? -0x00400000L : 0x00400000L);
> }
>
> if (y<0){ abs_y=-y; }
>
> if (x>=0){
> r = fixdiv(x - abs_y, x + abs_y);
> angle = FIX90 - fixmul(FIX90,r);
> }else{
> r = fixdiv(x + abs_y, abs_y - x);
> angle = FIX270 - fixmul(FIX90,r);
> }
> if (y < 0) return(-angle);
> else return(angle);
> }
>
> After some testing of the above, I get:
> fixatan2(50,50) = 90; // Should be 45
> fixatan2(80,20) = 144; // Should be 75.96
> fixatan2(20,80) = 36; // Should be 14.03
> fixatan2(-50,50) = -90; // Should be -45
>
> Can anybody see what I have done incorrectly? BTW, compiling on GCC for
> an ARM9 processor. No HW floats :) - thus all fixed.
>
> Thanks in advance,
> Peter
>
>
>
> This message was sent using the Comp.DSP web interface on
> www.DSPRelated.com
atan2(u,v) returns angles from 0 to (but not including) 360 degrees,
whereas atan(u/v) returns angles from -90 to +90 degrees. With either
function, it is necessary to compute an arctangent. I don't see one in
your code. One way to implement atan2() is getting atan(abs(u).abs(v)),
then translating the result to the proper quadrant determined by sgn(u)
and sgn(v). Where in your code to you get the numeric value of atan()?
Jerry
--
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������
Reply by pea●April 15, 20052005-04-15
Hello all,
Sorry to drag up old posts, but I found the original link and have been
trying to implement atan2 in fixed point. AFter it not working, I found
this post about the same algorithm and thought it a good place to start.
I'm using 16.16 fixed point. My version of the atan2 from the original
link follows. Please note - FIX90 and FIX270 are defines for PI/4 and
3*PI/4 in degrees. Looking at the code I thought it didn't matter that I
use degrees rather than radians, as I am expecting the result in degrees
anyway. Am I correct?
fixed fixatan2(fixed y, fixed x){
fixed abs_y = y;
fixed r, angle;
if (x==0){
if (y==0) return 0L;
else return ((y < 0) ? -0x00400000L : 0x00400000L);
}
if (y<0){ abs_y=-y; }
if (x>=0){
r = fixdiv(x - abs_y, x + abs_y);
angle = FIX90 - fixmul(FIX90,r);
}else{
r = fixdiv(x + abs_y, abs_y - x);
angle = FIX270 - fixmul(FIX90,r);
}
if (y < 0) return(-angle);
else return(angle);
}
After some testing of the above, I get:
fixatan2(50,50) = 90; // Should be 45
fixatan2(80,20) = 144; // Should be 75.96
fixatan2(20,80) = 36; // Should be 14.03
fixatan2(-50,50) = -90; // Should be -45
Can anybody see what I have done incorrectly? BTW, compiling on GCC for
an ARM9 processor. No HW floats :) - thus all fixed.
Thanks in advance,
Peter
This message was sent using the Comp.DSP web interface on
www.DSPRelated.com
Reply by ●February 1, 20052005-02-01
If you have specific questions let me know - I am the original author
of the atan2 algorithm on dspguru.
And the guys here are correct, the C code is easily converted to
fixed-pt. In fact it was first implemented in fixed-pt assembly, but
publishing DSP assembly code probably wouldnt get the point across to a
wide audience - hence the santized C version.
Reply by ●February 1, 20052005-02-01
If you have specific questions let me know - I am the original author
of the atan2 algorithm on dspguru.
And the guys here are correct, the C code is easily converted to
fixed-pt. In fact it was first implemented in fixed-pt assembly, but
publishing DSP assembly code probably wouldnt get the point across to a
wide audience - hence the santized C version.
Reply by JS●February 1, 20052005-02-01
If you have specific questions let me know - I am the original author
of the atan2 algorithm on dspguru.
And the guys here are correct, the C code is easily converted to
fixed-pt. In fact it was first implemented in fixed-pt assembly, but
publishing DSP assembly code probably wouldnt get the point across to a
wide audience - hence the santized C version.