Renormalizing the Sine Table
Phil Pilgrim (PhiPi)
Posts: 23,514
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:
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
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:
[b]PUB[/b] cos(x)
'' Cosine of the angle x: 0 to 360 degrees == $0000 to $2000,
'' with result renormalized to $1_0000.
[b]return[/b] sin(x + $800)
[b]PUB[/b] sin(x) : value | t
'' Sine of the angle x: 0 to 360 degrees == $0000 to $2000,
'' with result renormalized to $1_0000.
[b]if[/b] (x & $fff == $800)
value := $1_0000
[b]else[/b]
[b]if[/b] (x & $800)
t := -x & $7ff
[b]else[/b]
t := x & $7ff
value := [b]word[/b][noparse][[/noparse]*SINE][noparse][[/noparse]*t] - (t > 707)
[b]if[/b] (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
Comments
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
My Prop Info&Apps: ·http://www.rayslogic.com/propeller/propeller.htm
That's not a bad idea! It would only take 64 longs after all.
-Phil
[b]PUB[/b] cos(x) '' Cosine of the angle x: 0 to 360 degrees == $0000 to $2000, '' with result renormalized to $1_0000. [b]return[/b] sin(x + $800) [b]PUB[/b] sin(x) : value | t '' Sine of the angle x: 0 to 360 degrees == $0000 to $2000, '' with result renormalized to $1_0000. [b]if[/b] (x & $fff == $800) value := $1_0000 [b]else[/b] [b]if[/b] (x & $800) t := -x & $7ff [b]else[/b] t := x & $7ff value := [b]word[/b][noparse][[/noparse]*$e000][noparse][[/noparse]*t] + ((sin_corr[noparse][[/noparse]*t >> 5] >> (t & $1f)) & 1) [b]if[/b] (x & $1000) value := -value [b]DAT[/b] sin_corr [b]long[/b] $00000000,$00000000,$00208010,$30002008,$4001c000,$00401040,$a0420004,$24040a80 [b]long[/b] $0e042110,$21038000,$00084922,$30c4094a,$88600fc0,$4a150a48,$0000e312,$5529331c [b]long[/b] $07c31129,$aa92461f,$c07866d2,$4aa96cc3,$8e000732,$73256d49,$5ad9c1f0,$9ffe3369 [b]long[/b] $f32d55b3,$7b5b3c00,$671ff1d9,$fff3b5ad,$3db6bdb9,$dbeb6780,$55b3c01e,$7dcff9db [b]long[/b] $ef81e6d5,$bffddbfe,$ffeedfdb,$ffddebfe,$ffb6fb6f,$feedf6df,$ff77ffdf,$ffb7bdbf [b]long[/b] $23ff77df,$c1f6fedf,$bfdffff7,$ffffffdf,$defdedff,$ffef7feb,$ffdfff7f,$ffffdeff [b]long[/b] $ffdfffff,$fffeffff,$7fffffef,$ffffffff,$ffefefff,$ffffffff,$fffffeff,$fffff7f7 [b]long[/b] $ffffffff,$ffffffff,$fffff7ff,$bfffffff,$ffffffff,$ffffffff,$ffffffff,$ffffffff-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·==·$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·&·$800)
····t·:=·-x·&·$fff······'t ranges·from·$800 to $001
··else
····t·:=·x·&·$7ff······ 't ranges from $000 to $7FF
··value·:=·word[noparse][[/noparse]$e000][noparse][[/noparse]t]·+·((sin_corr[noparse][[/noparse]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······$07c31129,$aa92461f,$c07866d2,$4aa96cc3,$8e000732,$73256d49,$5ad9c1f0,$9ffe3369
··············long······$f32d55b3,$7b5b3c00,$671ff1d9,$fff3b5ad,$3db6bdb9,$dbeb6780,$55b3c01e,$7dcff9db
··············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,1 'Nothing to see here, Folks.
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Chip Gracey
Parallax, Inc.
-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
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.
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
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
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
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
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.
-Phil
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
My Prop Info&Apps: ·http://www.rayslogic.com/propeller/propeller.htm
$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
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