Trig. problem for an MCU

Your ratio a:b looks like the tangent of the angle, so perhaps a series expansion of the arctangent function might do.

Here's one that's pretty good for angles up to +/- 45 degrees:

f(x) = x - (1/3)x^3 + (1/5)x^5 - (1/7)x^7

equal,

cool thanks, that looks like it might be what I need, just a few multiplies too,

I spent a while looking at it but didnt see that, sometimes non simple trig and such stuff makes my brain want to be somewhere else seems to hapen more often these days.

Colin =^.^=

That formula, from the Maclaurin series, doesn't converge fast. Take a close look at the error at 45 degrees (x=1). The value is then:

This is 41.47 degrees, which is 3.53 degrees off-beam. You mentioned a need for "resolution of a degree." This won't cut it.

So expand it further. How many terms??

Well, let's look at a simple loop formulation of the above equation first, so we can test this:

Let fx = 0 For i = #_of_terms to 1 step -1 Let fx = 1/(2*i-1) - x^2*fx Next i Let fx = x * fx

For "1 - 1/3 + 1/5 - 1/7", we set #_of_terms to 4. From this, and looking at the worst case which is 45 degrees, we can now just plug in and see.

Turns out, you need 15 terms just to get to within .95 degrees of an angle of PI/2 and 29 terms to get within .49 degrees. Nasty. It doesn't converge well, at all. You might project from the above that it would take 58 terms to get within .25 degrees and you'd be right. It's that bad.

Continued fractions, which I think was used as insight by Euler for a different forumlation (which I won't place here -- look it up here:

if you want to see it) can be used, instead.

The following formulation is not present on the mathworld link above and I suspect it's not widely examined on the web. It does use division, but gets you there pretty fast so it may be worthwhile to consider it. You can make your own choices regarding what cycles out best on the dsPIC. But I'll present it for your examination:

x / (1 + x^2 / (3 + 0.8 * x^2))

For the x=1 case, this is 1/(1+1/3.8), which is about .36 degrees high.

I've picked the worst case to compare, though. 45 degrees is a problem for the Maclaurin series. If you can live with a somewhat smaller range of angles than 0-45, it turns out that those 4 terms in the Maclaurin series can get you within 1/2 degree at 38 degrees. (At

Hope this helps some.

Jon

P.S. I don't gather from your earlier comments that seeing this derived would help you much, so I won't bother writing out the details of the general form of the continued fraction for ATAN. Suffice it that the above formulation does have a derivation from basic rules.

However, the continued fraction version with one less step is:

x / (1 + x^2 / 3)

and with one more step is:

x / (1 + x^2 / (3 + 4 * x^2 / (5 + (9/7) * x^2)))

If you have a continued fraction like this, the rule for reducing from a larger number of steps to fewer, is to just treat the last x^2 term as 0 and then simplify. So the above expanded version would then be handled:

x / (1 + x^2 / (3 + 4 * x^2 / (5 + (9/7) * x^2))) x / (1 + x^2 / (3 + 4 * x^2 / (5 + (9/7) * 0))) x / (1 + x^2 / (3 + 4 * x^2 / (5 + 0))) x / (1 + x^2 / (3 + 4 * x^2 / 5)) x / (1 + x^2 / (3 + (4/5) * x^2)) x / (1 + x^2 / (3 + 0.8 * x^2))

Which is the formula I gave.

In addition to series solutions, this is the description of the Fortran function ATAN2D(a,b); if you have a suitable library of math functions, it might be just waiting there.

And since there are only 46 possibilities, you COULD just make a table (tand (0), tand (0.5),tand(1.5), ... tand(45.5)) and search the table for the largest value not greater than your ratio. The searching of a table is simpler than a long list of IF statements...

If divide is slow, consider that

a/b > c

is equivalent to

a> b*c

for any positive 'b'. Binary search of a table of 46 values involves 6 or 7 tests, so you might want to consider a no-divide solution.

b,

somewhere

Many thanks for the in depth reply, absolute accuracy isnt that essential as long as there arnt too many kinks such as where the quadrant 1/8s join up, basicaly its more increased resolution im after to avoid rotor jumps at low speed, the simplest 3 term is probably still beter than a short table lookup anyway. the dspic does actually have multiply accumulate, wich ive not made use of yet.

My input signals arent exactly that accurate anyway, but thats another problem, and im not often lucky enough for the errors to cancel each other out lol. at higher speeds the angle is determined ok by rate and time, but it would be usefull to just use 1 method.

the 3.5 deg error at 45' looks like it might make a bit of a kink, otherwise I would probably ignore it, but I could probably tweak the constants so the quadrants line up smoothly at the cost of accuracy elsewhere.

the sin=3/4-x^2+x approximation is used to convert the (shifted) angle back again, im not to sure of the actual acuracy but there seem to be no kinks in it.

the above sin approx was only one of several possibles I found each one having different flavours of error, with minimums of either offset, first second or third harmonics etc, I chose this as it was the quickest.

im not sure how it was derived as it doesnt look like traditional trig series im (vaguely) familiar with. I was wondering if there was a similar alternative to the usual series. I think it was derived using more numerical techniques.

I wasnt quite sure what I should be looking for before, a google for 'arctan approximation' threw up a few interesting things wich I will wade through, found a few graphs of errors for various polynomials etc.

thanks Colin =^.^=

Yes although the maplab gcc library does have atan its very long, but im not sure it has an atan2d equivalent. I already avoided the divides the way you suggest, and a binary search of a lookup table was my next step before I thought id stop and look at a more arithmetic aproach, as I found such an easy sin aproximation I thought it worth while looking, although as you point out 45 entries isnt very big.

thanks Colin =^.^=

I had forgotten to also add that since this Maclaurin series has an alternating sign, you can also consider averaging the alternating too-high and too-low values that result from adding more terms. In the above case, you could modify the snipped Maclaurin series by considering the last term in it as over-doing a good thing and instead taking just 1/2 of it (the last term) as a 'good measure,' so as to read:

f(x) = x - (1/3)x^3 + (1/5)x^5 - (1/14)x^7

This will be much better than before and takes the exact same time to calculate.

Jon

Aha so I might only need to work out what the last term multiplier has to be to get exactly 45' at x=1. that sounds promising.

Im surprised someone hasnt published modified optimised constants for short series, to give optimised results for lowest types of errors, harmonics etc.

a search showed up a couple of floating point implmentations of arctan2d in java group wich il have a closer look at . Here's a quick approximation, accurate to about .05 radians (in C): double arctan2( double y, double x ) { const double ONEQTR_PI = M_PI / 4.0; const double THRQTR_PI = 3.0 * M_PI / 4.0; double r, angle; double abs_y = fabs(y) + 1e-10f; if ( x < 0.0f ) { r = (x + abs_y) / (abs_y - x); angle = THRQTR_PI; } else { r = (x - abs_y) / (x + abs_y); angle = ONEQTR_PI; } angle += (0.1963f * r * r - 0.9817f) * r; if ( y < 0.0f ) return( -angle ); else return( angle );

Colin =^.^=

Or, alternatively, just use the arctangent.

Mark

I doubt that his MCU has a built-in atan function, and even if it did it would probably be too slow.

For a microcomputer and assuming anough flash/rom space, the fastest solution is probably a lookup table. Using a 16 bit sine table with values scaled ^ 2 so that the largest value just fits, then you use shifts, rather than */ to recover the final value. You only need to store one quadrant to cover both sine and cosine and a 16 bit table needs only 180 bytes for 1 degree steps.

Any other method is just hard work...

Chhris

final code :- (using 1/12 in last term)

// full circle = 16 bits

int atanii_(int a,int b) //a and b in range +-0x100 { int x=(a1); //last 7 bits are fraction int xx=x*x>>7; //x=+-0x80 int r=x*82; x=(x*xx)>>7; r-=x*27; x=(x*xx)>>7; r+=x*16; x=(x*xx)>>7; r-=x*7; return r; //r=+-0x2000 }

int atani(int a,int b) { int x; if(b>=a) if(b>=-a) x=atanii_(a,b); else x=0xc000-atanii_(b,a); else if(a>=-b) x=0x4000-atanii_(b,a); else x=0x8000+atanii_(a,b); return x; }

seems to progress smoothly accros semi quadrants estimated to take a couple of microseconds to execute. thanks to all who contributed.

now all I need to do is work out how to magnetise my rotor so it gives a reasonable sinewave on the hall sensors, or to better position the sensors.

Colin =^.^=

Oops. Sorry, my brain is taking a holiday this weekend ...

Mark