search

otopp / libfixmath

Automatically exported from code.google.com/p/libfixmath
0 stars 0 forks source link

Atan inaccurate #3

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago 
I am having problems with the atan implementation.

float x = .5;
float atanx = fix16_to_float(fix16_atan(fix16_from_dbl(x)));
printf("printf atan(x) = %.3f\n",atanx);

OUTPUT:
printf atan(x) = 0.000

It appears that since x has the decimal value of 32768 which is greater than 
29736 it goes into the if statement of _fix16_atan().

assuming inValue = .5 (32768 in fixed point)
tempOut = (1 + fix16_mul(inValue, inValue));    tempOut = 1 + .5*.5 = 1.25
tempOut = (1 + fix16_sqrt(tempOut));                 tempOut = 1 + sqrt(1.25) = 
2.118 (138805 in fixed point)
tempOut = (inValue / tempOut);                            tempOut = 
32768/138805 = 0
tempOut = _fix16_atan(tempOut);                        tempOut = _atan(0) = 0

So it appears that the division line breaks it.

Assuming I pass a number that in decimal is less than 29736
float x = .25;
float atanx = fix16_to_float(fix16_atan(fix16_from_dbl(x)));
printf("val = %ld\n",fix16_from_dbl(x));
OUTPUT:
printf atan(x) = 0.255

So even with numbers less than decimal 29736 the output is still wrong

It appears that the approximation
y = x + (x^3)/3 + (x^5)/5 + (x^7)/7 + (x^9)/9 + (x^11)/11
Is a much better approximation of tangent than the atan function from graphing 
it.

Original issue reported on code.google.com by joe.scha...@gmail.com on 22 Feb 2011 at 5:28

GoogleCodeExporter commented 9 years ago 
Very well spotted, I'm going to put that error down to my retarded state of 
mind late at night and the fact that I've never tested this. It seems that I 
forgot I was doing fixed point and mixed in some integer maths there, the 
correct code for those lines is:
tempOut = (fix16_one + fix16_mul(inValue, inValue));
tempOut = (fix16_one + fix16_sqrt(tempOut));
tempOut = fix16_div(inValue, tempOut);
tempOut = _fix16_atan(tempOut);
return (tempOut << 1);

The reason that if statement is there is because the approximation is quite 
wrong for numbers over that value. As for the inaccuracy of  the numbers below 
that threshold, I may need to lower the threshold or improve the accuracy of 
the algorithm. I did work out what accuracy should be needed but my maths was 
clearly very sketchy at the point when I wrote this specific function.

I'll commit it when I get home and get a chance to which may be later this 
week, in the mean time you seem quite interested in this project and have a 
good level of understanding, so your welcome to committer access if you want to 
join the project (don't worry I won't expect anything off you but it just means 
that you can fix any bugs that you may find).

Original comment by Flatmush@googlemail.com on 22 Feb 2011 at 5:26

GoogleCodeExporter commented 9 years ago 
I've just noticed a similar set of problems in asin, this should be:
fix16_t fix16_asin(fix16_t inValue) {
        if((inValue > fix16_one) || (inValue < -fix16_one))
                return 0; // Should cause error case.
        fix16_t tempOut;
        tempOut = (fix16_one - fix16_mul(inValue, inValue));
        tempOut = fix16_div(inValue, fix16_sqrt(tempOut));
        tempOut = fix16_atan(tempOut);
        return tempOut;
}

Original comment by Flatmush@googlemail.com on 22 Feb 2011 at 5:37

GoogleCodeExporter commented 9 years ago 
Looking at it, it might make sense for me to do the atan function in 64-bit 
(32.32) format to maintain accuracy unless I can find a faster way to get 
accurate results.

Original comment by Flatmush@googlemail.com on 22 Feb 2011 at 5:41

GoogleCodeExporter commented 9 years ago 
Ok, so I tried  out a  approximation I found at
http://www.dspguru.com/dsp/tricks/fixed-point-atan2-with-self-normalization

I have attached both implementations, the simple one which is <.08 radians off 
~ 4 degrees

the complex one which is < .01 radians off ~ .057 degrees 

They need some optimization obviously but the idea is there. 

We could use half angle approximations to find asin and acos?
http://en.wikipedia.org/wiki/Inverse_trigonometric_functions

let me know what you think.

Original comment by joe.scha...@gmail.com on 23 Feb 2011 at 12:44

Attachments:

GoogleCodeExporter commented 9 years ago 
So I did some tests on my ARM Cortex-M0

Attached is my code

With neither FLOAT_TEST or FIXED_TEST defined code size is 1741 Bytes
With FLOAT_TEST defined code size is 13176 Bytes
With FIXED_TEST defined code size is 3764 Bytes ( this is with the 4KB array 
removed, and the associated if)
With FLOAT_TEST and FIXED_TEST defined code size 13624 Bytes

to add the floating point sin it takes 11500 Bytes
to add the fixed point sin it takes ~2000 bytes
so a ~9500 Bytes savings.

to compute the floating point sine it took "Conversion Took 9210 ticks"
to compute the fixed point sine it took  "Conversion Took 2562"

Original comment by joe.scha...@gmail.com on 23 Feb 2011 at 2:36

Attachments:

GoogleCodeExporter commented 9 years ago 
Nice work, we should put these results on the wiki.

I hadn't thought of using this in such an embedded environment, I put in the 4K 
(16KiB) caches for speed because I was using this for games/3d on a less 
embedded environment, but it's clear that we should have some macros to disable 
them on smaller environments.

As for your code we should look into committing it, a few quick things to note:
You should work out the constants using other constants which can be reduced at 
compile time or using a calculator rather than using fix16_from_dbl, for 
example:
3 * (pi / 4) = (fix16_pi * 3) >> 2
or
3 * (pi / 4) = 3 * (3.1415926535897932384626433832795 * 65536) / 4 = 154415

Original comment by Flatmush@googlemail.com on 23 Feb 2011 at 9:31

GoogleCodeExporter commented 9 years ago 
Ok, I've committed the above changes to the code tree.

As for the complex atan2, I've made some changes to your code, tell me what you 
think. But I consider this to be your work and your contribution so since you 
now have committer access you should add this change in since it's you who 
deserves the credit for it.

A few things that I changed:

1. Use fix16_sdiv (saturated divide) so you don't have to worry about divide by 
zero, sdiv will simply return fix16_max if you do a divide by zero so it works 
well with most functions.

2. Like I said above, the code will be more accurate and faster if you use the 
proper constants for pi, so I replaced all the constants.

3. Made the code a bit more concise and inkeeping with the current code, though 
I kept your variable names the same since they're much nicer than the naming 
scheme I originally used. I'm no longer such a fan of putting temp infront of 
every local variable.

Original comment by Flatmush@googlemail.com on 23 Feb 2011 at 10:09

Attachments:

GoogleCodeExporter commented 9 years ago

Original comment by Flatmush@googlemail.com on 23 Feb 2011 at 10:10

GoogleCodeExporter commented 9 years ago 
Ok, my code had an error (the >> 16 in a constant).

I've run your test code with your version of the code and the version with a 
few small updates, and I can confirm that the updated version is very slightly 
more accurate. Included is the test with both versions of the code for 
comparison, the updated version is now fixed.

Original comment by Flatmush@googlemail.com on 23 Feb 2011 at 10:23

Attachments:

GoogleCodeExporter commented 9 years ago 
Hi,
Can you please upload the latest libfixmath after all the above fixes. 

BR,

Original comment by birdla...@gmail.com on 23 Feb 2011 at 12:49

GoogleCodeExporter commented 9 years ago 
Latest revision is uploaded.

Original comment by Flatmush@googlemail.com on 23 Feb 2011 at 2:31

GoogleCodeExporter commented 9 years ago 
It seems that the uploaded version does not contain the correct version on 
atan2. Can you plz check if I am not wrong.

Original comment by birdla...@gmail.com on 23 Feb 2011 at 4:55

GoogleCodeExporter commented 9 years ago 
I was going to let joe commit that change since it was his algorithm, but since 
you want it now I've uploaded a version which I've implemented in fixed point 
with no external calls, it's tested across the full range of fixed point 
numbers for accuracy though there may be some room for me to optimize it by 
removing some of the rounding that is required due to the dynamic range.

Original comment by Flatmush@googlemail.com on 23 Feb 2011 at 5:15

GoogleCodeExporter commented 9 years ago 
No problem with commiting I you seem to have improved it too.  It also seems 
that we are in far different timezones...

I will get a chance to test it in a little bit.  However It looks good, I like 
the no external calls part of it, also the low usage of division operations.  
The only part I personally dont like is the 64bit numbers, mainly because I 
don't know how well GCC handles 64bit numbers on a 32bit system, I assume it 
would handle them very well.  
But I will run a test when I get home using my "complex_atan2.c" and your 
updated atan2, and check codesize, conversion time differences, and size 
comparisons to the math.h atan2.  

Flatmush if you want to start the wiki, I could add an embedded section to it.  

Also is their a mailing list so I dont have to reply to the issue every time?

Original comment by joe.scha...@gmail.com on 23 Feb 2011 at 8:09

GoogleCodeExporter commented 9 years ago 
Not yet, I only just updated the project and it's not my only project. I 
actually didn't expect that many people to be interested in it really so that's 
a nice surprise. I guess I should look at sorting out a mailing list but I'm a 
bit rubbish at these admin things which is why I like googlecode.

Yeah the 64-bit integer arithmetic could be a problem, but without doing 
specific assembler for each platform (could do with gcc macros) there's no way 
of making it platform independent.

The real thing is that despite doing a lot of embedded stuff, this maths code 
has never really been the bottleneck of any of my projects, with graphics stuff 
it's almost always the triangle filler, etc. I'll have a look into benchmarking 
and optimizing it but if you beat me to the punch with that then be sure to 
commit it.

My timezone is GMT (UTC + 0), I'm guessing you're in America?

Original comment by Flatmush@googlemail.com on 23 Feb 2011 at 8:15

GoogleCodeExporter commented 9 years ago 
Mailing list created.

Original comment by Flatmush@googlemail.com on 23 Feb 2011 at 8:28

GoogleCodeExporter commented 9 years ago 
hi 

The latest version is more accurate. It is working fine with gcc but only one 
problem is that one of the compilers I use does not support int64. Is there any 
way around. 

Any benchmarking results for the library?

Original comment by birdla...@gmail.com on 24 Feb 2011 at 9:57

GoogleCodeExporter commented 9 years ago 
I'm going to make an issue ticket for 64-bit support, it's possible we can make 
some speed gains with it at the same time.

Original comment by Flatmush@googlemail.com on 24 Feb 2011 at 12:52