PDA

View Full Version : Renormalizing the Sine Table

Phil Pilgrim (PhiPi)
09-15-2009, 03:20 AM
The values in the Propeller's internal ROM sine table are normalized to \$ffff. In other words normalized(1) == \$ffff, and anything less than 1 is multiplied by \$ffff and rounded to the nearest integer to get the normalized value.

This isn't the way I would have done it. To me it makes more sense for \$1_0000 to equal one. That way sines and cosines can be used with fixed point multiplication and yield the expected results. Granted, this would mean that some values in the table would be too big to fit in a word, but that would have been okay. Just set them to zero and use a comparison on the table index to know when to set bit 16 of the looked-up value.

But that's water over the dam; the question is how to fix it. The table values are under their "true" \$1_0000-normalized values by one, 1300 times out of 2048. If we could roughly characterize those table indices for which the values are too low and add one to the result just for those, we could reduce the error. As a rule, the values at the upper end of the table are low more often than those in the lower part. So it makes sense to find an index threshold above which we add one to the result. I've done so and have found that a threshold of 707 on the index minimizes the number of errors to 176 too low by one and 216 too high by one, for a total of 392, or 30% of the original number.

Here's the Spin code that implements this adjustment:

PUB cos(x)

'' Cosine of the angle x: 0 to 360 degrees == \$0000 to \$2000,
'' with result renormalized to \$1_0000.

return sin(x + \$800)

PUB sin(x) : value | t

'' Sine of the angle x: 0 to 360 degrees == \$0000 to \$2000,
'' with result renormalized to \$1_0000.

if (x & \$fff == \$800)
value := \$1_0000
else
if (x & \$800)
t := -x & \$7ff
else
t := x & \$7ff
value := word[*SINE][*t] - (t > 707)
if (x & \$1000)
value := -value

It would be interesting to see if further refinements can be made, for example, by considering neighboring table values to determine which direction the value in question was rounded. But this simple correction, at least, is an improvement.

-Phil

localroger
09-15-2009, 03:42 AM
The error is never more than 1 part in 65536 or 0.001%, in the circle that is about 1/200 of a degree or 18 arc-seconds; this is several times finer than the best-case human visual acuity. It's very unlikely you would ever notice the difference in any RL application, and doing it this way allowed the Prop to spare you the nuisance of figuring out when 0 is 0 and when it's 1.

Rayman
09-15-2009, 04:29 AM
I suppose you could just make a 1-bit per element table for every value and thereby get the correct result...

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
My Prop Info&Apps: ·http://www.rayslogic.com/propeller/propeller.htm

Phil Pilgrim (PhiPi)
09-15-2009, 05:06 AM
Rayman,

That's not a bad idea! It would only take 64 longs after all.

-Phil

Phil Pilgrim (PhiPi)
09-15-2009, 05:36 AM
Thanks to Rayman's suggestion, the fix is in! Here's the code that does a 100% renormalization:

PUB cos(x)

'' Cosine of the angle x: 0 to 360 degrees == \$0000 to \$2000,
'' with result renormalized to \$1_0000.

return sin(x + \$800)

PUB sin(x) : value | t

'' Sine of the angle x: 0 to 360 degrees == \$0000 to \$2000,
'' with result renormalized to \$1_0000.

if (x & \$fff == \$800)
value := \$1_0000
else
if (x & \$800)
t := -x & \$7ff
else
t := x & \$7ff
value := word[*\$e000][*t] + ((sin_corr[*t >> 5] >> (t & \$1f)) & 1)
if (x & \$1000)
value := -value

DAT

sin_corr long \$00000000,\$00000000,\$00208010,\$30002008,\$4001c000, \$00401040,\$a0420004,\$24040a80
long \$0e042110,\$21038000,\$00084922,\$30c4094a,\$88600fc0, \$4a150a48,\$0000e312,\$5529331c
long \$ef81e6d5,\$bffddbfe,\$ffeedfdb,\$ffddebfe,\$ffb6fb6f, \$feedf6df,\$ff77ffdf,\$ffb7bdbf
long \$23ff77df,\$c1f6fedf,\$bfdffff7,\$ffffffdf,\$defdedff, \$ffef7feb,\$ffdfff7f,\$ffffdeff
long \$ffdfffff,\$fffeffff,\$7fffffef,\$ffffffff,\$ffefefff, \$ffffffff,\$fffffeff,\$fffff7f7
long \$ffffffff,\$ffffffff,\$fffff7ff,\$bfffffff,\$ffffffff, \$ffffffff,\$ffffffff,\$ffffffff

-Phil

cgracey
09-15-2009, 07:07 AM
Phil,

I think you might be able to get rid of the initial 'IF' structure by ANDing with \$FFF, instead of \$7FF. You'd need to add one more '1' bit·to your table.

PUB·cos(x)

''·Cosine·of·the·angle·x:·0·to·360·degrees·==·\$000 0·to·\$2000,
''···with·result·renormalized·to·\$1_0000.

··return·sin(x·+·\$800)

PUB·sin(x)·:·value·|·t

''·Sine·of·the·angle·x:·0·to·360·degrees·==·\$0000· to·\$2000,
''···with·result·renormalized·to·\$1_0000.

··if·(x·&·\$800)
····t·:=·-x·&·\$fff······'t ranges·from·\$800 to \$001
··else
····t·:=·x·&·\$7ff······ 't ranges from \$000 to \$7FF
··value·:=·word[­\$e000][­t]·+·((sin_corr[­t·>>·5]·>>·(t·&·\$1f))·&·1)
··if·(x·&·\$1000)
····value·:=·-value

DAT

sin_corr······long······\$00000000,\$00000000,\$00208 010,\$30002008,\$4001c000,\$00401040,\$a0420004,\$24040 a80
··············long······\$0e042110,\$21038000,\$00084 922,\$30c4094a,\$88600fc0,\$4a150a48,\$0000e312,\$55293 31c
··············long······\$ef81e6d5,\$bffddbfe,\$ffeed fdb,\$ffddebfe,\$ffb6fb6f,\$feedf6df,\$ff77ffdf,\$ffb7b dbf
··············long······\$23ff77df,\$c1f6fedf,\$bfdff ff7,\$ffffffdf,\$defdedff,\$ffef7feb,\$ffdfff7f,\$ffffd eff
··············long······\$ffdfffff,\$fffeffff,\$7ffff fef,\$ffffffff,\$ffefefff,\$ffffffff,\$fffffeff,\$fffff 7f7
··············long······\$ffffffff,\$ffffffff,\$fffff 7ff,\$bfffffff,\$ffffffff,\$ffffffff,\$ffffffff,\$fffff fff,1 'Nothing to see here, Folks.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔

Chip Gracey
Parallax, Inc.

Phil Pilgrim (PhiPi)
09-15-2009, 09:14 AM
Thanks, Chip. I think you're right if x is confined to the specified range. But I wanted to encompass negative angles and those over 360 degrees, too. In such cases, anding with \$fff would result in an index too big for either the ROM table or the correction table. However, I do need to double-check that out-of-range xes still return the correct results.

-Phil

Addendum: I've now tested the "\$7ff" version of the program with negative and out-of-range parameters, and it works okay. -P.

Post Edited (Phil Pilgrim (PhiPi)) : 9/15/2009 2:43:51 AM GMT

cgracey
09-15-2009, 02:16 PM
Phil,

consider this code:

··if·(x·&·\$800)
····t·:=·-x·&·\$fff······'t ranges·from·\$800 to \$001
··else
····t·:=·x·&·\$7ff······ 't ranges from \$000 to \$7FF

If (x & \$800) is true, and you·AND·(-x) with \$FFF, you'll get values ranging from \$800 (the extreme case) down to \$001. \$800 is the last valid table index, with the data word appearing at address \$F000..\$F001. The sine table is·sequentially read·from \$000 to \$7FF, then from \$800 back down to \$001, with that pattern repeating for negated lookup values to make the 2nd 180 degrees. This funny lookup pattern insures symmetrical zero-crossings and maximas. There should be no problem with negative angles coming in - those upper bits are always clipped off, anyway.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔

Chip Gracey
Parallax, Inc.

Phil Pilgrim (PhiPi)
09-15-2009, 03:01 PM
Chip,

Ah, so! You're right! The only time bit 11 is set for both x and -x is when x & \$fff == \$800.

I'm curious, though, since you've been very gracious in responding to my thread: when planning the sine table, was there ever a dilemma in your mind between normalizing to \$ffff vs. \$1_0000? I can see obvious advantages both ways. Assuming you had considered both, what was it that tilted your decision in favor of \$ffff?

Thanks,
-Phil

Post Edited (Phil Pilgrim (PhiPi)) : 9/15/2009 8:07:12 AM GMT

Rayman
09-15-2009, 05:51 PM
Phil: Can you also just multiply by 65536 and divide by 65535?

I think this error is pretty trivial until you start multiplying the results together, like in 3-d rotations or something...

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
My Prop Info&Apps: ·http://www.rayslogic.com/propeller/propeller.htm

Phil Pilgrim (PhiPi)
09-15-2009, 11:25 PM
Rayman,

No, it doesn't work because the numbers in the table are rounded from their true, fractional values. I tried a simulation of that in Perl and still got a lot of errors. Your correction table idea, along with Chip's simplification, are the gold standard here.

-Phil

cgracey
09-16-2009, 05:30 AM
Phil Pilgrim (PhiPi) said...

...when planning the sine table, was there ever a dilemma in your mind between normalizing to \$ffff vs. \$1_0000? I can see obvious advantages both ways. Assuming you had considered both, what was it that tilted your decision in favor of \$ffff?

Phil,

I just figured that in order·to get the most out of those words in the lookup table, I should push them up against their limits (maximum value = \$FFFF).

When I thought of normalizing the table to a whole binary number, I assumed I would have to settle on \$8000 for +/-1.0, effectively sacrificing·one·bit of·lookup-value resolution.

Even if there was no bit loss in having a 2^n limit (like +/-\$10000), the beauty of having a 2^n-1 limit (like +/-\$FFFF) is that the sign bit·can be flipped and with·a single SHR instruction, you've got a value that can be output directly to a·binary·DAC. You'd be in a mess if you had an 8-bit DAC (symmetrical off-center swing limited to +/-\$7F) and a value of +/-\$80 to output - you'd have to perform some scaling reduction. Also, a 2^n-1 limit affords maximum use of registers and memory, as there is never a possibility of extrema requiring one more bit of storage, effectively forcing all your computations to be one bit less in quality, in anticipation of the rare occasion when you'll have that leading '1' with many '0's trailing it.

It never occurred to me to do·a correction table like you made. But, it's probably because I'm biased towards DSP apps where those sine values are always getting scaled down by some other ugly numbers, and·~17 bits of resolution is already over the limit of what can likely be gleaned by playing analog games over the Propellers digital I/O pins.

For purists, though, I can see the interest in a correction·to +/-\$10000. This might be of value in some systems where only internal numbers are being handled. In practice, I don't·think it would make any discernable performance difference to an application which crosses into the analog domain. I'm sure it could help some programmers sleep better at night, at least until·they realize that there are much greater sources of error creeping into·their app from elsewhere that cannot be controlled nearly as well as their sine function.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔

Chip Gracey
Parallax, Inc.

Post Edited (Chip Gracey (Parallax)) : 9/15/2009 10:35:30 PM GMT

Phil Pilgrim (PhiPi)
09-16-2009, 06:04 AM
Chip,

Thanks for taking the time to explain your decision process! I was focused rather more myopically on internal computations involving rotations and Goertzel coefficients, both of which involve fixed-point multiplication. But, considering the points you've mentioned, along with the (now) relative ease of renormalizaiton when it's necessary, I can see that the strategy you chose does make the most sense. http://forums.parallax.com/images/smilies/smile.gif

-Phil

Rayman
09-18-2009, 07:21 PM
Phil: Can one just add a "1" to the sine table and then shift right by 1 to get your table with one less bit?

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
My Prop Info&Apps: ·http://www.rayslogic.com/propeller/propeller.htm

Tracy Allen
09-18-2009, 11:26 PM
The sine table is already mapping a 11 bit domain into a 16 bit range, so won't errors due to interpolation be much larger error than that one bit? For example, suppose you need the value PI/3 in the calculation. The value taken from the table without interpolation will not approximate that well. Or, suppose the result of one sin2 calculation feeds into the next angle, as in a sequence of rotations, aren't the values are usually going to hit between the 2 11 values so that the error will accumulate?

\$ffff2 = \$FFFE0001
which is quite different from \$100000000, but not so bad when reduced back to 16 or 11 bits. It could be a problem in a sum of squares with accumulation of error.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com (http://www.emesystems.com)

Phil Pilgrim (PhiPi)
09-19-2009, 12:19 AM
Tracy,

In my stereo 3D driver demo, I learned the hard way that you can't accumulate small transformations using the sine table without introducing some pretty "interesting" distortions. So, instead, I keep track of the net transformation, then compute new coordinates from the original, untransformed object after each small change.

I hadn't thought about it until you brought it up, but it would be interesting to write a sine function using a higher-resolution domain along with interpolation from the ROM table. The points are probably close enough that linear interpolation would yield adequate results, although even cubic spline interpolation would not be that difficult.

-Phil