DSPRelated.com
Forums

Trick: atan2() accurate linear interpolation

Started by Simon Hosie March 21, 2004
I have failed to find any evidence of this technique, but my research
skills are absolutely pathetic.  I don't suppose it's terribly likely
that I've managed to come up with a unique method, but on the off chance
that I have I'd rather see it in the public domain than anywhere else...
with my name on, if that would help me find a job in the field.

It's more than twice the speed of my Athlon's FPU while still being
accurate to at least 1/16384 of a circle.. I suppose that's good.


THIS WORK IS PLACED IN THE PUBLIC DOMAIN

Name: fixed-point atan2(); sub-quadrant folding for improved accuracy

Category: Algorithmic

Application: phase angle calculations

Advantages: fast/accurate compromise

Introduction:

Beyond folding the cordinates into a single quadrant so that they could
be passed to atan() without losing the direction of the vector, they can
trivially be folded even further to make linear approximation more
accurate.

The Trick:

At each conditional section in the following pseudo-C, we consider
whether the vector is outside of the half that we would like to reduce
it to, and if it is then we rotate the vector by a constant angle, and
offset the result accordingly.  Initially we do this by rearranging
variables and adjusting signs, but beyond 90 degrees it becomes
mathematically simpler to multiply by a vector (2^n, -1).  These vectors
produce non-ideal angles, but trivial arithmetic.

Each successive angle is slightly more than half the one before, which
is enough to ensure that one rotation will bring any possible value
within the new range, while the tangent of the angle may be expressed as
a binary shift.


	int i_atan2(int y, int x)
	{
		int result = 0;
		int y2;

		if ((x | y) == 0)
			return 0;	/* give up */

		if (y < 0)	/* if we point downward */
		{
			result += -180_deg;
			y = -y;
			x = -x;
		}
		if (x < 0)	/* if we point left */
		{
			result += 90_deg;
			y2 = y;
			y = -x;
			x = y2;
		}
		if (y > x)	/* 45 degrees or beyond */
		{
			result += 45_deg;
			y2 = y;
			y -= x;
			x += y2;
		}
		if (2 * y > x)	/* 26.565 degrees */
		{
			result += 26.565_deg;
			y2 = y;
			y = 2 * y - x;
			x = 2 * x + y2;
		}
		if (4 * y > x)	/* 14.036 degrees */
		{
			result += 14.036_deg;
			y2 = y;
			y = 4 * y - x;
			x = 4 * x + y2;
		}
		if (8 * y > x)	/* 7.125 degrees */
		{
			result += 7.125_deg;
			y2 = y;
			y = 8 * y - x;
			x = 8 * x + y2;
		}

		/* linear interpolation of the remaining 64-ant */
		result += (7.125_deg * 8) * y / x;
		return result;
	}


Example with 65536 units per circle, leaving optimisation to the
compiler.  The maximum error should be less than 2 units.  Accuracy can
be reduced or increased by changing the preprocessor directives.  Note
that there will be some error from the rounding of the result
adjustment, and so adding stages is not sufficient to improve accuracy
beyond a certain point.  Also note that the vector can grow in length
very quickly, and the input range shouldn't be limited. 

	int i_atan2(int y, int x)
	{
		int result = 0;
		int y2;

		if ((x | y) == 0) return 0;
		if (y < 0) { result += -32768; y = -y, x = -x; }
		if (x < 0) { result += 16384; y2 = y; y = -x; x = y2; }
		if (y > x) { result += 8192; y2 = y; y -= x; x += y2; }
		if (2 * y > x)
		{
			result += 4836;	/* atan(1/2) * 32768/M_PI */
			y2 = y;
			y = 2 * y - x;
			x = 2 * x + y2;
		}
	#if 1
		if (4 * y > x)
		{
			result += 2555;	/* atan(1/4) * 32768/M_PI */
			y2 = y;
			y = 4 * y - x;
			x = 4 * x + y2;
		}
	#if 1
		if (8 * y > x)
		{
			result += 1297;	/* atan(1/8) * 32768/M_PI */
			y2 = y;
			y = 8 * y - x;
			x = 8 * x + y2;
		}
	#if 1
		if (16 * y > x)
		{
			result += 651;	/* atan(1/16) * 32768/M_PI */
			y2 = y;
			y = 16 * y - x;
			x = 16 * x + y2;
		}
	#if 0
		if (32 * y > x)
		{
			result += 326;	/* atan(1/32) * 32768/M_PI */
			y2 = y;
			y = 32 * y - x;
			x = 32 * x + y2;
		}
		return result + 10427 * y / x;	/* up to atan(1/32) */
	#else
		return result + 10417 * y / x;	/* up to atan(1/16) */
	#endif
	#else
		return result + 10377 * y / x;	/* up to atan(1/8) */
	#endif
	#else
		return result + 10221 * y / x;	/* up to atan(1/4) */
	#endif
	#else
		return result + 9672 * y / x;	/* up to atan(1/2) */
	#endif
	}

Simon Hosie wrote:
> [...] Also note that the vector can grow in length very quickly, and > the input range shouldn't be limited.
"the input range *should* be limited". Post comes before proofread in the dictionary. Why not online?
"Simon Hosie" <simon.hosie@clear.REM-urrtaueght-OVE.net.nz> wrote in message
news:5622j1-6fu.ln1@bovine.muck.net.nz...
> [...] > I don't suppose it's terribly likely > that I've managed to come up with a unique method, but on the off chance > that I have I'd rather see it in the public domain than anywhere else... > with my name on, if that would help me find a job in the field. > [...]
Good for you! The method you describe is well known, and has been significantly refined. The magic word to Google for is CORDIC. Don't let that get you down -- it's a clever trick regardless of how many people discover it.
Matt Timmermans wrote:
> The method you describe is well known, and has been significantly > refined. The magic word to Google for is CORDIC.
Aww poop. Am I free to use it (in these modern times...), or will they confiscate my brain for even thinking of it?
Simon Hosie wrote:

> Matt Timmermans wrote: > >>The method you describe is well known, and has been significantly >>refined. The magic word to Google for is CORDIC. > > > Aww poop. Am I free to use it (in these modern times...), or will they > confiscate my brain for even thinking of it?
My HP-35 pocket calculator, which cost $400, used CORDIC to do trigonometry. I was told that the algorithm was implemented mechanically in the WW II Norden bombsight, but my guess is that it was used in another device and the statement results from confusion. Whatever, it's been around long enough for you to use it freely. The best description and explanation of CORDIC I know is a paper on Ray Andraka's web site, http://www.andraka.com/files/crdcsrvy.pdf. Search for CORDIC also at http://www.dspguru.com. Jerry -- Engineering is the art of making what you want from things you can get. &#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;
In article qps2j1-p0v.ln1@bovine.muck.net.nz, Simon Hosie at
sh1usenet@clear.REM-hararurhea-OVE.net.nz wrote on 03/22/2004 05:46:

> Matt Timmermans wrote: >> The method you describe is well known, and has been significantly >> refined. The magic word to Google for is CORDIC. > > Aww poop. Am I free to use it (in these modern times...), or will they > confiscate my brain for even thinking of it?
:-) it reminds me of something that Bob Adams used to say about whenever you think you invented something new, better check the Bell System Journals from the 50s or 60s. welcome to the club, Simon. (you have to pay Bill Gates $0.0000000000001 for every binary bit that you use.) r b-j
"Simon Hosie" <sh1usenet@clear.REM-hararurhea-OVE.net.nz> wrote in message
news:qps2j1-p0v.ln1@bovine.muck.net.nz...
> Matt Timmermans wrote: > > The method you describe is well known, and has been significantly > > refined. The magic word to Google for is CORDIC. > > Aww poop. Am I free to use it (in these modern times...), or will they > confiscate my brain for even thinking of it?
Reinventing CORDIC is still pretty impressive. Don't feel too bad about it. -Kevin
Jerry,

From you that means a lot.  Thanks!  Cordic was first published by Volder in
1956 I think.  Don't know if the WWII bombsights used a similar but
mechanical algorithm or not.  If so, it may not have been binary.  More
likely it used gears to work up sine and cosines similar to the ballistics
computers on WWII battleships.

Jerry Avins wrote:

> > The best description and explanation of CORDIC I know is a paper on Ray > Andraka's web site, http://www.andraka.com/files/crdcsrvy.pdf. Search > for CORDIC also at http://www.dspguru.com. > &#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;
-- --Ray Andraka, P.E. President, the Andraka Consulting Group, Inc. 401/884-7930 Fax 401/884-7950 email ray@andraka.com http://www.andraka.com "They that give up essential liberty to obtain a little temporary safety deserve neither liberty nor safety." -Benjamin Franklin, 1759
Ray Andraka wrote:

> Jerry, > > ... Don't know if the WWII bombsights used a similar but > mechanical algorithm or not. If so, it may not have been binary. More > likely it used gears to work up sine and cosines similar to the ballistics > computers on WWII battleships. > > Jerry Avins wrote:
That sounds much more reasonable than what I was told and mindlessly parroted. Ball-and-disc integrators were high tech then. Jerry -- Engineering is the art of making what you want from things you can get. &#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;