PDA

View Full Version : Some CORDIC documentation



Graham Stabler
10-10-2006, 04:22 AM
I need to do some real time trig for a differential odometer thing I'm working on, as such I've been struggling to get to grips with the CORDIC method. I finally got it and so for my own benefit I added these comments to the Memsic object that uses the Cordic method to compute angles given x and y lengths.

I have read quite a few explanations of CORDIC and found it hard going, in my explanation I hope to explain the process before filling in the maths and going on the realization. I've probably not done a better job than anyone else but it has helped me.

I'd welcome input from beginners and experienced folk on the accuracy and clarity of what I have written. In particular I'm unsure about the bits with setting the sign in the assembly section, note the ????????????.

Cheers,

Graham



' Documented Cordic
'
' A vector from the origin to a point may be represented as X and Y co-ordinates. And those co-ordinates
' if the vector is of unity length are equal to the sin and cosines of the angle between the x axis and the point.
'
' x = cos(Theta)
' y = sin(Theta)
'
' So we can see that if we find values of x and y corresponding to a given theta we have sin and cos(theta)
' and if we find the change in theta required to reduce y to zero we have found the angle theta. The CORDIC
' algorithm does just that.
'
'************************************************* ************************************************** *******************
'
' --- Performing a sin(theta) or Cos(theta), general description
'
' What we have are two vectors, one is the vector we want to determine the angle of or the x and y
' co-ordinates of (to get sine and cos) and the other is our estimate vector with co-ordinates x'and y'.
'
' In the case of finding Cos or Sin(theta) we start with our estimate vector on the x axis, it is then rotated by 45
' degrees towards the vector to be measured. 45 degrees is added to the angle accumulator and the values for
' x' and y' are modified. How do you do that? How can we work out where the end of the estimate vector is just by the
' angle when we can't do sin/cos? The trick is to choose the angles very carefully so we know what the change in x'
' and y' will be, more on that in the maths. The important word is change, we have to go thorough a process of tweaking
' the values of x' and y', estimate by estimate until they are correct.
'
' Say for instance that the angle we are trying to get the sin/cos of is +30 degrees, our rotation of the
' estimate vector has overshot. So the next rotation must be in the opposite direction. How big is the rotation?
' well it is smaller than before, 26.56 degrees, seems random but its not as we will see. So after this rotation
' the angle accumulator has 18.44 degrees in it and the values of x' and y' have been updated too, next the angle
' to be rotated is 14 degrees and the rotation direction is positive again. Now the angle accumulator
' has 18.44+14 = 32.44 degrees, getting warmer and the values in x' and y' are getting closer to sin(theta) and
' cos(theta) too.
'
' This process continues for a fixed number of iterations, the angle of rotation gets smaller each time as the correct
' angle is homed in on. At the end of the process you have your values of x' and y' and hence your sin(theta)
' and cos(theta) values.
'
'
'************************************************* ************************************************** *******************
'
' --- Getting the angle - General description.
'
' You can do a very similar thing to find the angle given the x and y co-ordinates and this is what the assembly
' code below does. To find an angle the values of x' and y' are set to those provided. The angle accumulator
' is set to zero and a similar process ensues. The difference in this case is that rather than aiming for a target
' vector the algorithm is trying to rotate the estimate vector to the x-axis, once it has done that the angle
' accumulator will contain the correct angle. As before it is performed rotation by rotation, the sign of y'
' determines the direction of rotation if it becomes positive the vector must be rotated CW and if negative CCW
' to make sure it moves towards the x-axis. So although not used directly the values of x' and y' are still
' important in deciding if to add or subtract the increasingly small rotation angles to the angle estimate.
'
'
'************************************************* ************************************************** *******************
'
' --- Some maths
'
' As we know a vector can be represented by a pair of co-ordinates x,y
'
' If that vector is rotated around the origin by an angle theta the new co-ordinates of its endpoint become
'
' x' = x cos(theta) - ySin(theta) ' From some trig (look up on net if interested)
' y' = y cos(theta) + xSin(theta)
'
' As Tan(theta) = Sin(theta)/Cos(theta) then it may be re written as
'
' x' = cos(theta)( x - yTan(theta))
' y' = cos(theta)( y + xTan(theta))
'
' This looks a little worse if anything, we now have a Tan multiplied by y. If only yTan(theta) could be easily
' calculated. Well in binary, multiplies/divides by a power of 2 (2,4,8...) can be done by bit shifting. To multiply by 2 you
' simply shift the bits 2^1 places to the left, 0010 => 0100 voila! Similarly you can divide with left bit shifts so if Tan(theta)
' was made a power of 2 the yTan(theta) could be done with simple shifts.
'
' As binary itself shows us you can represent a number by adding up lots of progressively smaller bits so this is what
' is done in the Cordic algorithm. We choose Tan(theta) to be a limited set of powers of 2, to be precise we set:
'
' Tan(theta) = 2^-i (note: 2^-i means 1/(2^i) and i is the iteration number,)
'
' when i = 1 Tan(theta) = 2^-1 and theta = 45 degrees
' when i = 2 Tan(theta) = 2^-2 and theta = 26.56 degrees
' ....
'
' Sound familiar?
'
' We can easily put these angles in a look up table and access them as we go along, the example below
' uses 20 of them, that's not many for some amazing accuracy.
'
' Doing some more mathematical substitution and defining Xi as the current estimate and Xi+1 as the next
' we have:
'
' Xi+1 = Cos(atan(2^-i)).[Xi - YiDi(2^-i)]
' Yi+1 = Cos(atan(2^-i)).[Yi + XiDi(2^-i)]
'
' Its not as bad as it looks. Firstly I have dropped the ' for clarity and secondly atan(2^-i) has been subbed for
' theta because that is the angle we shall be rotating by. There is also a new variable called Di, I suppose that D
' might stand for decision as this value is determined by the direction of rotation required and may be +/-1.
'
' The whole thing can be neatened up by defining:
'
' Ki = Cos(atan(2^-i))
'
' Xi+1 = Ki.[Xi - YiDi(2^-i)]
' Yi+1 = Ki.[Yi + XiDi(2^-i)]
'
' Here's the clever part, Ki does not contain anything that varies other than i, and i always does the same thing,
' it goes from 1 : n. The thing to understand is that without Ki the estimate vector gets longer and longer and
' will end up being longer than one, that will mean that Xi and Yi will no longer represent Sine and Cosine of
' Theta. Ki simply scales the vector to unity and it does the same thing every time you run the cordic
'
' The final effect of Ki is that the result has been multiplied by K1*K2*K3*K4 ...... Kn
'
' This may be given the name K (the cordic gain) and applied to Xi and Yi eiher after or before iterating
' __________________________________________________ _________________________________________________
'| |
'| Aside: K, the Cordic gain |
'| |
'| K = Product(Ki) from K1-Kn ' This means K1*K2*Kn... |
'| |
'| K = Product(Cos(atan(2^-1))) ' Subbing for Ki |
'| |
'| K = Product(1/sqrt(1+2^-2i)) ' according to some boffin |
'| |
'| and that becomes 0.607252935 very quickly as n gets bigger (because cos of a small angle is one). |
' __________________________________________________ __________________________________________________ |

' So cutting a long story short you either multiply X and Y by 0.607252935 at the beginning or at the end if you
' care about there real lengths.

' So now we have the more dainty equations of:
'
' Xi+1 = [Xi - YiDi(2^-i)]
' Yi+1 = [Yi + XiDi(2^-i)]
'
' A reminder on how to use them.
'
' To find sin/cos:
'
' 1. Begin with Xi = 1 and Yi = 0 so the estimate is on the x-axis
' 2. Iterate the equation, updating the values of Xi and Yi, Increment/decrement the angle accumulator according to Di
' with the values in the look up table, 45, 26.5 etc. Work out the next value for D by comparing the the angle accumulator
' with the angle provided (is the estimate less or greater than the angle).
' 3. Repeat 2 for a given number of steps, more steps means greater accuracy.
' 4. Scale Xi and Yi by 0.607252935
' 5. And relax
'
' To find theta given x and y:
'
' 1. Begin with Xi = x and Yi = y , the values given.
' 2. Iterate the equation, updating the values of Xi and Yi. Increment/decrement the angle accumulator with the values in the
' look up table, 45, 26.56 etc decide on Di by looking at sign of Yi in order to rotate the estimate towards the x-axis.
' 3. Repeat 2 for a given number of steps
' 4. Relax, the angle accumulator holds the answer no need for scaling.
'
'
'************************************************* ************************************************** ***************
' And now some assembly code:
'
'
' This is code taken from the mesmic2125 object, it computes angle given x any y values. I have added my comments.
' Comments between lines tend to be explanatory of the process, those on the lines are about the code and may help
' those new to assembly (like me)

DAT 'cx and cy are the supplied variables.
'ca will hold the eventual answer


cordic ' The initial vector could be in any of the 4 quadrants, the basic algorithm will only work for angles between
' +/- 90 degrees but you want to be able to get an answer in any instance. You also want to ensure that the final
' answer is between 360 degrees.
'
' The number system for angle is carefully chosen, it varies from 0-360 or rather it varies from 0-0, as it overflows
' back to 0 when full. This means that if your starting vector is in the 4th quadrant (down and right) and is say
' 10 degrees from the x-axis if the angle is set to zero degrees at the start then as the estimate approaches the
' x-axis the angle becomes 340.
' __________________________________________________ __________________________________________________ ____________
'| |
'| Why? Its 2s-complement, |
'| |
'| take an 8 bit number 00000000 then take one from it you get 11111110 which if you forget about the fact it |
'| should be neagitive is 254 or more importantly 255-1 |
'| |
'| As a further note, in assembly the program has to keep track of whether a number is supposed to be signed |
'| or not. e.g. the first line of this code checks for negativeness (wc)so we must assume cx is being used as a |
'| signed variable. |
'|________________________________________________ __________________________________________________ ______________|

' If the vector is in the 1st quadrant the angle is just added as the estimate vector moves towards the x-axis.
'
' If in the 3rd quadrant (left and down) the angle is set to 180 at the start and the vector flipped to the 1st quadrant.
' This means that the final angle is 180 + the swept angle as it should be.
'
' Finally for the 2nd quadrant the angle is set to 180 at the start and the vector flipped to the 4th quadrant the final result
' becomes 180-swept angle.
'

' Make cx positve as the algoritm always operates in quadrants 1 and 2.
' This makes Q2=>Q1 and Q3=>Q4
abs cx,cx wc ' take the absolute of cx
' If was in Q3 or Q2 flip about x-axis
' So Q2=>Q4 and Q3=>Q1 overall.
if_c neg cy,cy ' negate cy if cx was negative

' Decide on starting angle, 0 or 180
mov ca,#0 ' load 0 into ca, the estimated angle
rcr ca,#1 ' Load with 180 if was on left

movs :lookup,#table ' Puts address of table into sumc command
mov t1,#0 ' Initialize t1 (angle shift 2^-t)
mov t2,#20 ' Initialize t2 (number of interations)

'The loop
:loop 'First calculate dX(i+1) => Yi.2^-t
mov dx,cy wc ' Load dx with cy cary is set if negative
sar dx,t1 ' Shifts dx by t1 while preserving sign
'Then calculate dY(i+1) => Xi.2^-t
mov dy,cx ' Load dy with cx
sar dy,t1 ' Shifts dy by t1 while preserving sign
'Now add to the estimate with sign dependant on sign of Yi
sumc cx,dx ' Add -dx to cx if c=1 or dx if c=0
sumnc cy,dy ' Add -dy to cx if c=0 or dy if c=1
'Next add or subract the specified rotation from the angle accumulator
:lookup sumc ca,table ' Add -table to ca if c=1 or table if c=0
' Remembering that the source field of the sumc command "table"
' is being over written to address different elements
' of the table below
add :lookup,#1 ' Increments the destination of lookup sumc command
add t1,#1 ' increment t1, this is i
djnz t2,#:loop ' decrement t2, keep looping if not zero.

cordic_ret ret ' return to rest of program (whatever it is)

'Table of atan(2^-i)


table long $20000000 'atan(1) = 45 degrees
long $12E4051E 'atan(1/2) = 26.56 degrees
long $09FB385B 'atan(1/4)
long $051111D4 'atan(1/8)
long $028B0D43 'atan(1/16)
long $0145D7E1 'atan(1/32)
long $00A2F61E 'atan(1/64)
long $00517C55 'atan(1/128)
long $0028BE53 'atan(1/256)
long $00145F2F 'atan(1/512)
long $000A2F98 'atan(1/1024)
long $000517CC 'atan(1/2048)
long $00028BE6 'atan(1/4096)
long $000145F3 'atan(1/8192)
long $0000A2FA 'atan(1/16384)
long $0000517D 'atan(1/32768)
long $000028BE 'atan(1/65536)
long $0000145F 'atan(1/131072)
long $00000A30 'atan(1/262144)
long $00000518 'atan(1/524288) = 1e-4 degrees

' Uninitialized data

t1 res 1
t2 res 1

dx res 1
dy res 1
cx res 1
cy res 1
ca res 1

Post Edited (Graham Stabler) : 10/10/2006 7:02:10 PM GMT

Beau Schwabe (Parallax)
10-10-2006, 05:21 AM
Graham Stabler (http://forums.parallax.com/member.php?u=46678),

Great job!

When I wrote the Memsic object I was at the Parallax office in Rocklin, and had Chip over my shoulder helping me
out.· Originally I had an object that used a reverse lookup table deriving Arcsine from of the existing Sine table in
ROM.· Although it worked it was not as fast, and didn't have the precision of Chip's CORDIC routine.· At the time
he happened to be using the CORDIC stuff for something else he was working on, so it was literally a cut-n-paste
from his object into mine with only a few minor tweaks here and there.· As for a detailed explanation on how the
CORDIC routine works, understand it a little.
·
·

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Beau Schwabe (mailto:bschwabe@parallax.com)

IC Layout Engineer
Parallax, Inc.

Tracy Allen
10-10-2006, 11:53 AM
Great job, Graham, seconded! Thanks for posting this.

I'll add a couple of comments.

1) This statement regarding the cosine term,

> At each stage of the process k does a known and consistant thing to the result.
> This allows K to be taken out of process and just applied at the end as a scaling factor.
> As it happens it very quickly converges to 0.607252935"

It is not k by itself, but the _product_ of terms K(i) that quickly converges:


{k(0) * k(1) * k(2) * ... * k(n-1) --> 0.607252935}


It converges quickly, because cosine approaches unity as the rotation angles get smaller and smaller.

Also, there is a central insight behind CORDIC is is fuzzy there. The rotation for each angle can be either + or -, and the cosine term also contains the factor, (D) for decision, but this is the insight...


k = Cos(Di * atan(2^-i)) <--- note the Di is included!

but from elementary trig...
cosine(theta)=cosine(-theta).

k = Cos(atan(2^-i)) <--- therefore the Di can be dropped



2) In the item, "To find theta, given x and y". It is true that there is no need for scaling of theta when the algorithm is finished. But there is a "furthermore". At the end, when the vector has been rotated to the x axis, the value cx now contains the magnitude of the original vector, which may be also be of interest. That needs to be scaled by the CORDIC gain.

3) It is really instructions like sumc and sumnc that set the Propeller apart from other processors. In one step, it accomplishes one of the core operations of CORDIC, to speedily add or subtract conditonal on a flag. My hat is off to Chip! He is always thinking ahead!

4) You asked about the business with the


rcr ca,#1 ' sets sign bit if cx was negative ??????????????


I don't quite understand it either, but I am sure it is very clever and there will be an aha!

Note that the value Chip has chosen to represent 45 degrees is $2000_0000. I say "chosen", because, other things being equal, that is indeed an arbitrary choice. You might want to choose decimal 45_000_000 to represent 45 degrees, and to construct the table, you would simply enter that starting number into excel and have it pop out the values of 45_000_000 * atan(2^-i) for i=0 to 31. It is simply a matter of scale and the CORDIC will work just fine at any scale. This table is constructed with the values, $2000_0000 * atan(2^-i). The full circle is 8 time 45 degrees, the full circle $ffff_ffff, or +/- $7fff_ffff. That of course gives the best possible angular resolution.

Now for my take on the question. The CORDIC algorithm has to be restricted to work with angles in the range of -90 to +90 degrees, or +/- $3fff_ffff. By priming the value of angle (ca) with the sign bit, the additions and subtractions of the algorithm are relected through the origin into the left half plane. These things always take a tight twist of logic! If the sign bit is set, then ca starts at $8000_0000, which is 180 degrees, or -180 degrees if you wish. Positive rotations add a positive value to that, say a 45 degrees moves to $A000_0000, which is 225 degrees in quadrant 3, or equivalently - 135 degrees. And a negative rotations move into quadrant 2.

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

cgracey
10-10-2006, 02:10 PM
Graham, Beau, and Tracy:

I've found that CORDIC scaling is best performed BEFORE the rotation so that the results never overflow the unit circle. This means multiplying X (and maybe Y) beforehand by the constant 0.607252935 that Tracy pointed out. This is, of course, a pain in the rear. In 32-bit quality, that number in binary is:

%0.1001101101110100_1110100110101000 ("_" separates groups of 16 bits)

For 16-bit quality, this would round up to:

%0.1001101101110101

This would take ~50 instruction cycles for a looped multiplier to compute.

There's a faster way, though. That sequence of 1's and 0's can be searched for sub-patterns which might have exploitable symmetries. I looked at that sequence and found a few sub-patterns which repeated, but %10001 was the most agreable since it occurs four times at equal relative offsets, allowing those four occurrences to be·summed in·two folded operations. The last occurrence requires individual attention since it·didn't fit into the first four and was at a different relative offset, anyway:



······· 0.1001101101110101·· <---·CORDIC voodoo·scaling term
······· ------------------
+ >>0···· 10001 ^^ ^^^ ^ ^·· <--- made up of goblins
+ >>3······· 10001 ^^^ ^ ^
+ >>6·········· 10001^ ^ ^
+ >>9 ············ 10001 ^
+ >>11·············· 10001·· <--- orphan



Here's some code that will scale·X by %0.1001101101110101 (0.607254) in just 12 instructions:


··············· sar···· x,#1··········· 'get %0.10001 of x
··············· mov···· t1,x··········· '(····· 1··· )
··············· sar···· t1,#4·········· '(+········ 1)
··············· add···· x,t1

··············· mov···· t1,x··········· 'get %0.10011001 of x
··············· sar···· t1,#3·········· '( ······· 10001)
··············· add···· t1,x··········· '(+·····10001·· )

··············· mov···· t2,t1·········· 'get %0.10011011011001 of x
··············· sar···· t2,#6·········· '( ·········· 10011001)
··············· add···· t2,t1·········· '(+·····10011001····· )

··············· sar···· x,#11·········· 'get %0.1001101101110101 of x
··············· add···· x,t2··········· '( ··············· 10001)
······································· '(+···· 10011011011001· )

·

{k(0) * k(1) * k(2) * ... * k(n-1) --> 0.607252935} must be one of those cosmically special numbers!

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


Chip Gracey
Parallax, Inc.

Post Edited (Chip Gracey (Parallax)) : 10/10/2006 8:15:29 AM GMT

Graham Stabler
10-10-2006, 03:58 PM
Tracey,

1) you are right about the k thing, it is K not k that converges, I'll correct that, bad wording.

2) As for the d in the cosine, if you define k as not having d in it i.e. you seperate the angle from its sign you get the same result and don't have to consider the symetry or lack of it. I've seen explanations on the web that mention the symmetry and those that don't. Not sure if what I have done is really wrong.

3) Interesting about the length of the vector being in Cx, I'll add that info.

4) There is nothing random about Chips choice for 45 degrees, if you multiply $2000_0000 by 8 which should give 360 degrees you get $1_0000_0000 which in 32bits is zero, and 360 and zero are the same thing. This provides I assume the maximum potential for resolution? Oh I think you said that.

5) I'm going to work out this Ca thing, I think you are probably right, I'll come back with more info.

Graham Stabler
10-10-2006, 04:05 PM
Chip,

1) Thanks, I think :) Now some new stuff to get my head around! I read it was more convenient to pre-scale so that the vector becomes of unit length at the end but as I didn't know how to do it, I went for the easy option in my explanation. I might add a part 2, a more verbose description of what you have explained and an example that calculates the cin and cosines of a given angle.

2) In a recent thread I asked about sources of clever tricks with binary, this is definitely one, can you give me some insight on where you learnt all this? And yes I know you were born with a microprocessor in your hand :)

Graham

Graham Stabler
10-10-2006, 07:14 PM
I have updated the description of K and have added more description to the assembly now I think I finally get it.

Chip, PLEASE check it!

Graham

Graham Stabler
10-10-2006, 07:20 PM
Chip/Tracey,

If you want to calculate the Sin of an angle you start with the estimate vector on the x axis and give it unit length, isn't multiplying 1 by 0.607252935 even easier :)

The only time you need to do the multiplication is if you are trying to work out the vector length and angle given x and y, so you are finding sqrt(x^2 + y^2)

is that correct?

Graham

Graham Stabler
10-10-2006, 08:05 PM
Corrected which quadrant it which.

cgracey
10-11-2006, 12:16 AM
Graham,

My incomplete view of the CORDIC algorithm is that it is good for two basic things:

1) Rotating an (X,Y) point around the origin by driving the desired angle towards 0.

····· a) Initializing X to Ro, Y to 0, and angle to Theta is useful for polar to cartesian conversion

···· ·b) Initializing X to 0.6072539, Y to 0, and setting angle is useful for COS/SIN generation - you get both at once!

····· c) Simply rotating (X,Y) by an angle is often useful for graphics and signal synthesis.

2) Clearing the angle and then rotating·(X,Y)·by driving Y towards 0 yields the original (X,Y) angle and length in X

····· a) This is effectively cartesian to polar conversion

Depending on·your needs, the algorithm can be simplified around the setup and takedown code which goes·before and after the rotation code. The·rotation code is always the same in function, though the the decision point·may be·either the sign of the angle or the sign of Y. You can make a loop or straight-line code it.

I've read that the CORDIC algorithm can also be used·for log conversions. I imagine the main difference would be replacing the lookup table with some other values. Tracy could probably figure out how to do this.

As far as numerical tricks go to simplify computations: they always revolve around the commutative and associative properties of addition, subtraction, multiplication, and division (which are the building blocks of integrals, differentials, and everything higher). Programming real-time signal-processing applications is so much fun because you're building a functional puzzle that must be sub-solved many times along the way to completion. It's always about putting more stuff in, and then finding ways to take·almost all of it back out, so that just a few instructions remain to perform the essential work.

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


Chip Gracey
Parallax, Inc.

Post Edited (Chip Gracey (Parallax)) : 10/10/2006 7:44:21 PM GMT

Bean
10-11-2006, 12:56 AM
Graham,
· Very nice. But I do have one question. In the text at the top you have:


' As Tan(theta) = Sin(theta)/Cos(theta) then it may be re written as
'
' x' = cos(theta)( x - yTan(theta))
' y' = cos(theta)( x - yTan(theta))


I can't see how X' and y' can be equal to the exact same formula ?

Bean.



▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Cheap used 4-digit LED display with driver IC·www.hc4led.com (http://www.hc4led.com)

Low power SD Data Logger www.sddatalogger.com (http://www.sddatalogger.com)
SX-Video Display Modules www.sxvm.com (http://www.sxvm.com)

Don't mistake experience for intelligence. And vis-vera.
·

Graham Stabler
10-11-2006, 02:10 AM
Bean, that was the special 45degree only version :)

Now corrected and thanks for taking the time to read!

Thanks Chip, I now have an almost complete version of the code that will provide Sin(theta) and Cos(theta) (at the same time yeh!), I'm just debugging it. Learning the binary tricks is good fun but it can be mind bending decoding other peoples clever code.

I think I may need that divide routine after all as I am doing some resolving, I think my vector length will need to be multiplied by 0.65.... So I can end up with r.cos(theta)

Cheers

Graham

Tracy Allen
10-11-2006, 02:45 AM
In the Stamp version of the CORDIC arctangent, www.emesys.com/BS2mathC.htm (http://www.emesys.com/BS2mathC.htm), I had to incorporate a wrapper to prescale numbers up to larger values. For example, if x=3 and y=3, the angle is should obviously be 45 degrees and the length sqr(6)=2.449489742783, but the algorithm would not converge well, unless the small x and y were both multiplied by a contant factor such as 2^10 * 0.607252935. You can't effectively prescale a small integer by a fraction like 0.607252935! By scaling both x and y by the same factor, the angle is unaltered, and the length scales by the same factor, and the scale factor can be compensated at the end. As I write this it's not obvious how the memsic Prop algorithm above wiil respond to small integers. Having 32 bits of resolution should converge very close to the correct angle, within
long $00000518 'atan(1/524288) = 1e-4 degrees
But I think a prescaling wrapper would be necessary for precision with "small" numbers on SQR(x^2 + y^2).


I didn't learn anything about CORDIC until I had an email exchange with Al Williams about his math coprocessor. (Where is Al these days? -- I haven't seen him on the forums for a long time). That was when I got interested in implementing (maybe _kluging_ would be a better description) it on the BASIC Stamp. Ray Adraka's article is the best I have seen for an overview of what can be accomplished with CORDIC and a succinct treatment of the procedures. The underlying machinery is the same. For exmple, to get to log and exponential, the table contains the values of the hyperbolic trig funcitons instead of the ordinary trig functions. As for creative use in signal processing, I think Chip has it down, bar none. Wow, that quick algorithm for multiplying by 0.607252935 is cool! I was wondering too how you came up with that decomposition.

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

Bean
10-11-2006, 03:59 AM
When I wrote an ATAN CORDIC program in SX/B, I scaled the inputs according to:
If either input (X or Y) was > 16384 then both were divided by 2.
Then...
If both inputs are < 4096 then both are multiplied by 2.

Here is the code if anyone is interested...

Bean.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Cheap used 4-digit LED display with driver IC·www.hc4led.com (http://www.hc4led.com)

Low power SD Data Logger www.sddatalogger.com (http://www.sddatalogger.com)
SX-Video Display Modules www.sxvm.com (http://www.sxvm.com)

Don't mistake experience for intelligence. And vis-vera.


Post Edited (Bean (Hitt Consulting)) : 10/10/2006 9:04:09 PM GMT

cgracey
10-11-2006, 07:43 AM
Tracy Allen said...
Wow, that quick algorithm for multiplying by 0.607252935 is cool! I was wondering too how you came up with that decomposition.
I just drew the pattern onto paper and started crossing out groups of bits. I was surprised to find that a couple of sub-patterns emerged right away. I chose %10001 because of the 4 occurrences at constant relative offsets. This makes me think that a searcher could be made which would scan a target multiplier pattern and come up with minimized solutions. Aside from adding, there's always subtracting, too. I suppose there are hundreds of solutions to some patterns. I wonder if any optimizing compilers ever bothered to do this kind of thing. Maybe even non-linear formulae could be realized, as well. This subject fascinates me. One the one hand you have some mathematical objective, while on the other you have a processor capable of doing only certain things. How to connect the two most efficiently is the challenge.

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


Chip Gracey
Parallax, Inc.

Post Edited (Chip Gracey (Parallax)) : 10/11/2006 2:58:02 AM GMT

Tracy Allen
10-11-2006, 07:58 AM
Here is an Excel spreadsheet that underscores how rapidly the cordic gain converges to 0.607252935.



I cos(atn(1/2^I)) product
0 0.70710678118655 0.70710678118655
1 0.89442719099992 0.63245553203368
2 0.97014250014533 0.61357199107790
3 0.99227787671367 0.60883391251775
4 0.99805257848289 0.60764825625617
5 0.99951207608708 0.60735177014130
6 0.99987795203470 0.60727764409353
7 0.99996948381879 0.60725911229889
8 0.99999237069278 0.60725447933256
9 0.99999809265682 0.60725332108988
10 0.99999952316318 0.60725303152913
11 0.99999988079073 0.60725295913895
12 0.99999997019768 0.60725294104140
13 0.99999999254942 0.60725293651701
14 0.99999999813736 0.60725293538591
15 0.99999999953434 0.60725293510314
16 0.99999999988359 0.60725293503245
17 0.99999999997090 0.60725293501477
18 0.99999999999272 0.60725293501035
19 0.99999999999818 0.60725293500925
20 0.99999999999955 0.60725293500897
21 0.99999999999989 0.60725293500890
22 0.99999999999997 0.60725293500889
23 0.99999999999999 0.60725293500888
24 1.00000000000000 0.60725293500888
25 1.00000000000000 0.60725293500888
26 1.00000000000000 0.60725293500888
27 1.00000000000000 0.60725293500888
28 1.00000000000000 0.60725293500888
29 1.00000000000000 0.60725293500888
30 1.00000000000000 0.60725293500888
31 1.00000000000000 0.60725293500888



The term cos(atn(1/2^n)) in the product can be re-written as follows, eliminating the cos and atn.

http://forums.parallax.com/attachment.php?attachmentid=43567

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

Post Edited (Tracy Allen) : 10/11/2006 1:07:45 AM GMT

Phil Pilgrim (PhiPi)
10-11-2006, 11:30 AM
Chip,

Automating the decomposition process intrigues me, too. As it turns out, it may not be that difficult — at least for 16-bit word sizes. The reason I say that is that the search space is fairly small, as combinatorial problems go, making an exhaustive search feasible.

Anyway, for a w-wide pattern of n ones, and w - n zeroes, there are comb(n, w) = n! / ((n-w)! w!) possible combinations. At each of the 17 - w shift positions this pattern can take, we have three choices: add, subtract, or do nothing. So there are 3 ** (17 - w) different ways to handle each pattern. Therefore, the total search space for a w wide pattern of n one bits is

····comb(n, w) × 317 - w

Summing this over a pattern width of n to 16 yields the following for each n from 2 to 8:




n | Search Size
2 | 48,427,344
3 | 24,212,652
4 | 12,102,756
5 | 6,042,096
6 | 3,002,484
7 | 1,472,070
8 | 699,570




For any given pattern and any number of adds/subtracts of that pattern, it's easy to compute both the resulting coefficient (to see if it's even correct) and, if it is, the number of machine instructions (the "cost") required to perform the computation. Then it becomes a cost-minimization problem over the entire problem space.

Okay, so checking 90 million different algorithms will probably take awhile. Whether or not to apply heuristics to decrease the search space will depend on how many times you have to do it, I guess. If I only had to do it once, I'd go for the exhaustive search. If that took too long and I had to do it more than once, I might look for a cleverer approach.

Unfortunately, the cost function is by no means continuous over the search domain. So gradient search methods are out. One could use genetic algorithms, but they don't guarantee the best solution, only a good one.

Maybe I'll try the exhaustive technique, just to see how long it takes.

-Phil

Post Edited (Phil Pilgrim (PhiPi)) : 10/11/2006 4:40:14 AM GMT

cgracey
10-11-2006, 01:21 PM
Phil Pilgrim (PhiPi) said...

Automating the decomposition process intrigues me, too. As it turns out, it may not be that difficult — at least for 16-bit word sizes. The reason I say that is that the search space is fairly small, as combinatorial problems go, making an exhaustive search feasible...

...For any given pattern and any number of adds/subtracts of that pattern, it's easy to compute both the resulting coefficient (to see if it's even correct) and, if it is, the number of machine instructions (the "cost") required to perform the computation. Then it becomes a cost-minimization problem over the entire problem space.

Phil,

This would work for multiplies and divides and it would be very handy. I'm anxious to hear all about anything you come up with.

I've daydreamed for a long time about some program that would help you write code for math functions, but not necessarily by knowing anything about math. It would·just simulate random instruction groupings in order to·arrive at·some desired characteristic. Say you want a random number generator with a certain distribution, or you need something to approximate a root only to a certain quality, or you want some fuzzy filter, etc. It would be like some kind of chemistry set for assembly coding. You could use it to synthesize odd little 5-instruction sequences or maybe even whole systems, depending on your behavioral criteria. I bet it would come up with stuff that would completely defy analysis, but (apparently) work fine. Maybe all that is practical is something that can deal with up to 8-instruction sequences, or so, but even that would be educational to play with. It could teach us some things.


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


Chip Gracey
Parallax, Inc.

Phil Pilgrim (PhiPi)
10-11-2006, 04:01 PM
Chip,

I remember reading about some guys in England who did that with an FPGA. They "trained" it, using a genetic algorithm, to act as a frequency discriminator. It worked brilliantly, and with a minimum of gates, even though the resulting gate assemblage made no intutive sense.

The "multiply decomposition" problem, and others like it, would be an interesting arena for trying a GA, since the "alleles" could be limited to a small subset of the available instructions, i.e. shifts, adds, subtracts, Booleans, etc. Ideally, the search space would be representated in such a way that small structural changes (resulting from crossover and mutation, say) would entail small behaviorial changes. This is almost never the case with digital computer instructions, though, where changing even a single bit can result in wildly different behavior. But, at least, eliminating the worst offenders (e.g. jumps and calls) can help to regularize the search domain and "smooth" the resulting utility (cost) function.

-Phil

James Long
10-11-2006, 09:41 PM
I hate to break in here with pure ignorance........I don't remember CORDIC in any of my past math.....although that was a long time ago.

What is CORDIC used for? A simple explaination would be great.·I'm just try to get this in the soft fleshy part contained in the skull.



James L

Phillip Y.
10-11-2006, 10:08 PM
I remember an article in the HP journal about the new HP35 calculator and that used CORDIC functions,
the function for TAN did not need to be scaled because the factor canceled out or reduced to below the 10 digit resolution,
instead of calculating each binary digit it was done in decimal one digit at a time.

I also reviewed an article buy someone who built a prototype CORDIC IC for the DOD (wo a multiplier),
the author mentioned that he used a trick he called "servo feedback" during the rotation calculations
to eliminate the error/scale factor, I think this amounted to combining/doing both calculations simultaneously.

The articles are buried in my archives.

Mike Green
10-11-2006, 10:10 PM
I can't speak about earlier uses (probably military), but HP used CORDIC algorithms in their ground-breaking programmable calculators.
These were serial BCD shift register processors with maybe 16 digits in a register, could add, subtract, and compare, had maybe 16 to 32
words of memory (16 digits each) originally with enough masked ROM for the algorithms. The CORDIC algorithms allowed them to provide
all of the trig and log functions to high precision.

Tracy Allen
10-11-2006, 10:18 PM
James,

If you want a good introduction to CORDIC, download the classic article from Ray Andraka's site. www.andraka.com/cordic.htm (http://www.andraka.com/cordic.htm).

I doubt if you learned it in school, because it is not usually taught, which to me seems very strange.

It is the basis for calculation of the trig and transcentental functions and inverses in practically every hand calculator. It is a very fast O(n) successive approximation algorithm that is well suited to small integer processors or FPGAs. It has its origins in navigation computers for avionics, where the trig functions and inverses have to be computed in real time. It has exotic applications in signal processing and filtering. I've heard speculation that mathematicians at the NSA have figured out ways to break encryption codes using CORDIC methods, but that would be a stretch.

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

Graham Stabler
10-11-2006, 10:43 PM
James, in case the above reference is a bit heavy for you CORDIC is used to work out Sin, cos, arcsin, arccos, tan, arctan and logs as well. It does this very quickly using a very small look up table and needs only bit shifts and additions.

The above algorithm takes an X and y co-ordinate and provides an angle accurate to about .0004 degrees in only 20 iterations and it doesn't have to do much in each iteration.

Gaham

Tracy Allen
10-11-2006, 10:48 PM
It's always worthwhile to study Chip's voodoo code, to see how he exorcises the goblins.

Here is a reinterpretation, a graphic decomposition of the cordic gain he posted yesterday.

http://forums.parallax.com/attachment.php?attachmentid=43577

The critical observation is that the factor 17=%10001 is repeated. (Phil and Chip, I'd be very interested too, to see how this process of decomposition, to discover those repeating factors, could be automated!?)

Chips' algorithm effectively applies parentheses around the groupings as in the bottom row, so it computes x * 17/(2^5), and then uses that grouping "a" to compute super grouping "b", and then uses supergrouping "b" another time, and finally adds in the orphan at the end. Here is Chip's pasm code to do that. Due to the folding, it takes 12 instructions and 12 clock cycles to multiply x times the prescaler.





sar x,#1 'get %0.10001 of x
mov t1,x '( 1 )
sar t1,#4 '(+ 1)
add x,t1

mov t1,x 'get %0.10011001 of x
sar t1,#3 '( 10001)
add t1,x '(+ 10001 )

mov t2,t1 'get %0.10011011011001 of x
sar t2,#6 '( 10011001)
add t2,t1 '(+ 10011001 )

sar x,#11 'get %0.1001101101110101 of x
add x,t2 '( 10001)
'(+ 10011011011001 )





Just to see how it would work out, I tried it with the parentheses and grouping as in the second line. That computes x*17/(2^5) as grouping a, but then repeatedly shifts and adds by 2^3, and then finally by 2^2 at the end. It uses only 11 instructions, but because of the loop, it takes 18 clock cycles instead of Chip's 12. (caveat, I'm not sure if my untried code is corrrect, and I'd appreciate comment if anyone spots an error.)




sar x,#1 'get %0.10001 of x
mov t1,x '( 1 )
sar t1,#4 '(+ 1)
add x,t1

mov t1,x ' x* 17/(2^5)
mov t2,#3 ' going to repeat 3 times

loop sar x,#3 ' add factor, x*17/(2^n), for n=8,11,14
add x,t1
djnz t2,#loop ' one clock cycle for loop back, 2 for fall through

sar x,#2
add t1,x ' add factor 17/(2^16)



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

Post Edited (Tracy Allen) : 10/11/2006 3:53:42 PM GMT

cgracey
10-11-2006, 11:39 PM
James Long said...

What is CORDIC used for? A simple explaination would be great.·I'm just try to get this in the soft fleshy part contained in the skull.


I remember reading that the CORDIC algorithm was developed in the 1950's to serve as·a practical in-flight navigation "computer".·The acronym stands for Coordinate Rotation DIgital Computer.
It makes me wonder how many other such ingenious algorithms may be possible, but haven't been discovered.

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


Chip Gracey
Parallax, Inc.

Harley
10-12-2006, 12:38 AM
CORDIC! Hadn't seen that term for decades.

Back in the mid-to-late '60s (yes, way past 'sell by date'), I worked on a military project for Missile Battery equipment. CORDIC was used by the discrete F/Fs and gates computer (BC = before computers or chips). Somewhere I may still have a document describing CORDIC. Seemed complex. Was used to correct for the non-flat earth when sites are distant from each other; their reference planes were tilted. Target data had to remove that 'aiming' error.

One of those 'been there, seen that' situations. Good to know CORDIC is still useful. http://forums.parallax.com/images/smilies/yeah.gif

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Harley Shanko
h.a.s. designn

Bean
10-12-2006, 01:48 AM
Here is a document that does a decent job of explaining CORDIC atan/hyp algorithm.

http://www.gmw.com/magnetic_measurements/Sentron/sensors/documents/AN_125KIT_manual.pdf#search=%22an125_kit%20cordic% 22

Bean.


▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Cheap used 4-digit LED display with driver IC·www.hc4led.com (http://www.hc4led.com)

Low power SD Data Logger www.sddatalogger.com (http://www.sddatalogger.com)
SX-Video Display Modules www.sxvm.com (http://www.sxvm.com)

Don't mistake experience for intelligence. And vis-vera.


Post Edited (Bean (Hitt Consulting)) : 10/11/2006 6:51:54 PM GMT

Graham Stabler
10-12-2006, 02:30 AM
Race you! http://forums.parallax.com/images/smilies/freaked.gif

Cordic FFT (http://www-dt.e-technik.uni-dortmund.de/publikationen/documents/Eusipco2004.pdf#search=%22cordic%20fft%22)

Tracy Allen
10-12-2006, 03:12 AM
This is one of those of those topics that needs to be revisited over and over again and from different approaches. Thank you for those references.

Here are some more:

www.andraka.com/cordic.htm (http://www.andraka.com/cordic.htm)
Ray Andraka Consulting
By far the most concise yet thorough discussion of the subject I know of. Download the pdf file, "Survey of CORDIC Algorithms for FPGAs".

www.dspguru.com/info/faqs/cordic2.htm (http://www.dspguru.com/info/faqs/cordic2.htm)
Grant R. Griffin dspGURU.
An informative FAQ, with additional links to the online and printed literature.


www.ee.ualberta.ca/courses/ee401/microboard/cordic_CCink.html (http://www.ee.ualberta.ca/courses/ee401/microboard/cordic_CCink.html)
Ingo Cyliax
CORDIC swiss army knife for computing math functions.. An article from Circuit Cellar INK. It even includes a program for calculating sine and cosine on the BASIC Stamp II. A different
approach from the one taken here.

www.ddj.com/documents/s=1070/ddj9010f/9010f.htm (http://www.ddj.com/documents/s=1070/ddj9010f/9010f.htm)
Pitts Jarvis
IMPLEMENTING CORDIC ALGORITHMS A single compact routine for computing transcendental functions. A tutorial from Dr. Dobbs Journal, October 1990.

people.csail.mit.edu/hqm/imode/fplib/cordic_code.html (http://people.csail.mit.edu/hqm/imode/fplib/cordic_code.html)
www.math.utep.edu/Faculty/helmut/presentations/1997/cordic/wcordic.html (http://www.math.utep.edu/Faculty/helmut/presentations/1997/cordic/wcordic.html)
Helmut Knaust
How Do Calculators Calculate? A tutorial from University of Texas.

Volder, J.E., 1959, The CORDIC Trigonometric Computing Technique,
IRE Transactions on Electronic Computers, V. EC-8, No. 3, pp. 330-334
seminal paper

www.math.ucl.ac.be/%7Emagnus/num1a/cordic.txt (http://www.math.ucl.ac.be/%7Emagnus/num1a/cordic.txt)
Vladimir BaykovInteresting newsgroup posting by Vladimir BAYKOV

www.cnmat.berkeley.edu/%7Enorbert/cordic/node3.html (http://www.cnmat.berkeley.edu/%7Enorbert/cordic/node3.html)
Norbert Lindlbauer, Berkeley Center for New MusicDiscussion of CORDIC architectures with a view toward music synthesizer

www.emesys.com/BS2mathC.htm (http://www.emesys.com/BS2mathC.htm)
Some tutorial on theory and code for sin/cos and for atn/hyp on the BASIC Stamp.
http://forums.parallax.com/showthread.php?p=563729
A discussion with Peter Verkaik about prescaling and about cordic for the Javalin Stamp.

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

Graham Stabler
10-12-2006, 03:36 AM
Tracy, I should have said earlier how much help your site was to me and I certainly have not finished reading!

I liked your spreadsheet showing the convergence of k so I went a step further and have produced a spreadsheet that will calculate SIN and COS of theta using CORDIC, I didn't do it out of philanthropy but because my damn program keeps cocking it up and I want to know when and why!

Find attached, the 360=>32bit column computes the decimal number the angle will appear if I output it to the TV (for debug), I need a similar thing for Xi and Yi possibly.

Graham

James Long
10-12-2006, 03:45 AM
Thanks everyone for the related links.....I'll have to do some reading.....interesting...we learn something new everyday.



James L

cgracey
10-12-2006, 04:07 AM
Graham Stabler said...
Race you! http://forums.parallax.com/images/smilies/freaked.gif

Cordic FFT (http://www-dt.e-technik.uni-dortmund.de/publikationen/documents/Eusipco2004.pdf#search=%22cordic%20fft%22)
What a great idea!

In computing an FFT you must multiply your samples by a sine and a cosine. This is the heart of the FFT and it takes 95% of the computing power - the rest is just adds and increments. CORDIC would simplify this central·operation into a single step. Just set X to your sample, clear Y, rotate by the angle, and you've got the two products in X and Y. The CORDIC expansion could be ignored, since it would be systemic.

CORDIC would do wonders for unraveling the·accummulations into angles and magnitudes, too. That would use the steer-Y-to-zero mode, instead of steer-angle-to-zero used to rotate X and Y.

In fact, this brings the FFT down to "understandable" (never mind the butterfly part), since you're dealing more in polar concepts which are usually obscured by such ill-fitting nomenclature as "real" and "imaginary". It took me forever to understand the root concept of·the·Fourier transform·because of how muddled the definition is due to how math is taught/explained. It's no wonder almost nobody can explain WHY it works, since their brain was prematurely clogged with the formality of the definition. Same goes for integrals and derivatives. I bet 99% of the people who teach calculus can't explain how come reducing an exponent by 1 and multiplying by the original exponent produces a derivative (I don't know, myself, but I'm really curious). Obviously, the person who figured it out understood it, with INTENT. Today, almost nobody knows WHY. They just know the recipe. Education ought to be about understanding WHY, not just HOW. It sure would make things easier for me, anyway. Fortunately, we got Paul Baker at Parallax now and he can "interpret" this math stuff pretty well. And Tracy Allen knows nearly everything, no matter how esoteric, or at least how to figure it out.

Here is one outstanding question that not even Tracy has answered (and it might be because there isn't an answer):

····· How can you sum polar coordinates (angle, distance) without converting back and forth·to/from cartesian? Anyone?

Anyway, Graham, this is a great idea!

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


Chip Gracey
Parallax, Inc.

Graham Stabler
10-12-2006, 04:15 AM
I don't know why I didn't do this before

Here is the new spreadsheet that will calculate sin(theta) and cos(theta) given theta and on the second sheet will provide Theta and vector length given X and Y

Find attached

Graham

Graham Stabler
10-12-2006, 04:24 AM
Chip, since someone mentioned FFTs a while back I have been thinking about it, that idea of the multifunction scope, logic analyser and spectrum analyser got me excited. Not to mention all the other ways you can use FFTs, a fastish photodiode and a laser pointer and you could probably measure blood flow (perfusion) by analyzing the spectrum for example.

When I finally got the cordic method I stuck "cordic fft" into google and got even more excited. I've been trying to refresh myself on FFT's, I used to use them a lot but its been a while.

You are dead right about the way things work being lost in the terminology, e has a lot to answer for.

Graham

cgracey
10-12-2006, 04:47 AM
Jumping the gun here·a bit, I'm thinking about something that is really needed badly:

A SIMPLE METHOD TO SUM POLAR VECTORS

Here's why: in signal processing we are always "slumming it" in cartesian land because of the simplicity of vector summation: you can simply add the X's and add the Y's. Meanwhile, we're having to lookup or generate sine and cosine, and scale each of them. A lot of work! When it comes time to rotate (X,Y), WHOAH! Thank goodness for CORDIC here. If we could reside in polar land, rotation would be just an add to the angle term. The length could be independently manipulated, too. Never any need to deal with a sine or cosine until we were generating the final output. But, because there's no apparent way to sum vectors in the polar domain, we are forced into the mud of cartesian processing.

If there was an efficient algorithm for summing polar vectors, in polar form, signal processing would become very simple to grasp. Signal processing needs polar treatment. Cartesian is an awful fit.·We need to "live" in polar land just to have the clarity of mind to think about signal processing.

I bet there's some undiscovered algorithm that would make this possible.


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


Chip Gracey
Parallax, Inc.

Post Edited (Chip Gracey (Parallax)) : 10/11/2006 10:13:54 PM GMT

cgracey
10-12-2006, 12:31 PM
Graham Stabler said...
Chip, since someone mentioned FFTs a while back I have been thinking about it, that idea of the multifunction scope, logic analyser and spectrum analyser got me excited. Not to mention all the other ways you can use FFTs, a fastish photodiode and a laser pointer and you could probably measure blood flow (perfusion) by analyzing the spectrum for example.

That's a fascinating idea. How does the spectrum reflect blood flow, though? It seems the beam is wide enough that it would shine through many blood cells which would probably en masse look like band-limited white noise surging with each pump of the heart. Is that how you would tell?

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


Chip Gracey
Parallax, Inc.

Tracy Allen
10-12-2006, 02:45 PM
Chip, regarding the polar addition of vectors. There might be a very constrained way to go about it. Here is a diagram that shows two vectors at angles A1 and A2, lengths L1 and L2,. You want to add them to find A3,L3.

http://forums.parallax.com/attachment.php?attachmentid=43590

The angle labeled a6 depends only on the angles A1 and A2 of the original vectors, not their lengths. However, the angles a7 and a8, and the resultant angle A3, all do depend on the lengths going in. If there is a constraint on the lengths, if for example they are both the same, or they are unit vectors, or they have some fixed ratio, then it is possible to solve it geometrically and the lengths drop out. For example, if L2 = L1, equal lengths, then,

A3 = A2 + (A1 / 2)

{correction: should be A3 = (A2 + A1) / 2}

Angle A3 depends on only A2 and A1, not L2 and L1 (except that the constraint has to be maintained as the vectors rotate). I haven't thought it through for all the quadrants.


The resultant length is still a problem. Observe again from the geometry that angle a6 depends only on A1 and A2, not on the lengths. if the lengths L1 and L2 are equal as above, then

L3 = 2 * L1 * sin (a6 / 2).

I'm sure that is not the general solution you are looking for. But it could be useful for something. I am thinking of the spirograph, where you see the pattern traced out by resultant sum of two vectors, where the ratio of the rotation frequecies and the ratio of the lengths are both constants.

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

Post Edited (Tracy Allen) : 10/12/2006 5:27:10 PM GMT

Peter Verkaik
10-12-2006, 03:17 PM
Tracy,
for equal lengths, shouldn't that be A3 = (A1+A2)/2
after all, adding (1,0) to (0,1) yields (1,1)
(0 + pi/2)/2 = pi/4

regards peter

Paul Baker
10-12-2006, 03:18 PM
Amusing Tracy, I've been working on the magnitude part. The magnitude of the addition of two polar vectors is sqrt(ra2 + rb2 + rarbcos(Θa-Θb)). The third term turns out to be the dot product of the two vectors.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Paul Baker (mailto:pbaker@parallax.com)
Propeller Applications Engineer
[/url][url=http://www.parallax.com] (http://www.parallax.com)
Parallax, Inc. (http://www.parallax.com)

Peter Verkaik
10-12-2006, 03:27 PM
Paul,
shouldn't that be
sqrt(ra2 + rb2 + 2rarbcos(Θa-Θb))
because adding 2 vectors of equals length and angle results
in a vector of double length with the same angle.

regards peter

cgracey
10-12-2006, 03:32 PM
Tracy Allen said...
Chip, regarding the polar addition of vectors. There might be a very constrained way to go about it. Here is a diagram that shows two vectors at angles A1 and A2, lengths L1 and L2,. You want to add them to find A3,L3.

http://forums.parallax.com/attachment.php?attachmentid=43590

The angle labeled a6 depends only on the angles A1 and A2 of the original vectors, not their lengths. However, the angles a7 and a8, and the resultant angle A3, all do depend on the lengths going in. If there is a constraint on the lengths, if for example they are both the same, or they are unit vectors, or they have some fixed ratio, then it is possible to solve it geometrically and the lengths drop out. For example, if L2 = L1, equal lengths, then,

A3 = A2 + (A1 / 2)

Angle A3 depends on only A2 and A1, not L2 and L1 (except that the constraint has to be maintained as the vectors rotate). I haven't thought it through for all the quadrants.


The resultant length is still a problem. Observe again from the geometry that angle a6 depends only on A1 and A2, not on the lengths. if the lengths L1 and L2 are equal as above, then

L3 = 2 * L1 * sin (a6 / 2).

I'm sure that is not the general solution you are looking for. But it could be useful for something. I am thinking of the spirograph, where you see the pattern traced out by resultant sum of two vectors, where the ratio of the rotation frequecies and the ratio of the lengths are both constants.



Well, maybe it's analogous to the first iteration of some CORDIC-like algorithm which could be an all-case solution. Thanks for thinking about this. I really believe that with enough thinking, we could come up with at least some algorithm, though not necessarily a processor-efficient one. I've searched the web for this solution, but I've seen on many sites "vectors cannot be added in polar form" (period). I pretty clearly remember having some complex formula for this, though, in my high-school physics class where we needed to add force vectors. If we did come up with·some formula we could start to think about how to break it down.

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


Chip Gracey
Parallax, Inc.

Paul Baker
10-12-2006, 03:56 PM
It's late, I'll show my math and see if you can figure if I missed something

rc=sqrt((racos(Θa) + rbcos(Θb))2 + (rasin(Θa) + rbsin(Θb))2)

= sqrt(ra2cos2(Θa) + rb2cos2(Θb) + rarbcos(Θa)cos(Θb) + ra2sin2(Θa) + rb2sin2(Θb) + rarbsin(Θa)sin(Θb))

= sqrt(ra2(cos2(Θa) + sin2(Θa)) + rb2(cos2(Θb) + sin2(Θb)) + rarb(cos(Θa)cos(Θb) + sin(Θa)sin(Θb)))

substituting cos2(Θ) + sin2(Θ) = 1

= sqrt(ra2 + rb2 + rarb(cos(Θa)cos(Θb) + sin(Θa)sin(Θb)))

substituting sin(Θa)sin(Θb) = (cos(Θa-Θb) - cos(Θa+Θb))/2

and cos(Θa)cos(Θb) = (cos(Θa-Θb)·+ cos(Θa+Θb))/2

= sqrt(ra2 + rb2 + rarb(cos(Θa-Θb)/2 - cos(Θa+Θb)/2 + cos(Θa-Θb)/2·+ cos(Θa+Θb)/2))

and with the final culmination of terms:

= sqrt(ra2 + rb2 + rarbcos(Θa-Θb))

Did I miss something somewhere?

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Paul Baker (mailto:pbaker@parallax.com)
Propeller Applications Engineer
[/url][url=http://www.parallax.com] (http://www.parallax.com)
Parallax, Inc. (http://www.parallax.com)

Post Edited (Paul Baker (Parallax)) : 10/12/2006 9:02:21 AM GMT

Peter Verkaik
10-12-2006, 04:00 PM
Chip,
how about this:

assume 2 vectors in the same quadrant
r1.angle(phi1) and r2.angle(phi2) and r1<r2

result = r1.angle(phi1) + r2.angle(phi2)
= r1.angle(phi1) + r1.angle(phi2) + (r2-r1).angle(phi2)
= r1*K1.angle((phi1+phi2)/2) + (r2-r1).angle(phi2)
= r1'.angle(phi1') + r2'.angle(phi2')
for the next iteration reassign r1,phi1,r2 and phi2 and repeat the above
until r2' reaches near zero (this determines the accuracy)

The difficulty is calculating r1'

regards peter

Post Edited (Peter Verkaik) : 10/12/2006 9:25:47 AM GMT

Peter Verkaik
10-12-2006, 04:04 PM
Paul, you missed the factor 2 in working out (a+b)*(a+b) = a*a + 2*a*b + b*b

It's late, I'll show my math and see if you can figure if I missed something

rc=sqrt((racos(Θa) + rbcos(Θb))2 + (rasin(Θa) + rbsin(Θb))2)

= sqrt(ra2cos2(Θa) + rb2cos2(Θb) + 2rarbcos(Θa)cos(Θb) + ra2sin2(Θa) + rb2sin2(Θb) + 2rarbsin(Θa)sin(Θb))

= sqrt(ra2(cos2(Θa) + sin2(Θa)) + rb2(cos2(Θb) + sin2(Θb)) + 2rarb(cos(Θa)cos(Θb) + sin(Θa)sin(Θb))

substituting cos2(Θ) + sin2(Θ) = 1

= sqrt(ra2 + rb2 + 2rarb(cos(Θa)cos(Θb) + sin(Θa)sin(Θb))

substituting sin(Θa)sin(Θb) = (cos(Θa-Θb) - cos(Θa+Θb))/2

and cos(Θa)cos(Θb) = (cos(Θa-Θb)·+ cos(Θa+Θb))/2

= sqrt(ra2 + rb2 + 2rarb(cos(Θa-Θb)/2 - cos(Θa+Θb)/2 + cos(Θa-Θb)/2·+ cos(Θa+Θb)/2)

and with the final culmination of terms:

= sqrt(ra2 + rb2 + 2rarbcos(Θa-Θb))

regards peter

Paul Baker
10-12-2006, 04:08 PM
Ok I see it now, thanks Peter.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Paul Baker (mailto:pbaker@parallax.com)
Propeller Applications Engineer
[/url][url=http://www.parallax.com] (http://www.parallax.com)
Parallax, Inc. (http://www.parallax.com)

Graham Stabler
10-12-2006, 04:18 PM
Chip Gracey (Parallax) said...


That's a fascinating idea. How does the spectrum reflect blood flow, though? It seems the beam is wide enough that it would shine through many blood cells which would probably en masse look like band-limited white noise surging with each pump of the heart. Is that how you would tell?


You get a doppler signal from the interference between coherent light comming from static tissue (unshifted) and moving tissue(the blood). As there is blood traveling at different speeds and in slightly different directions your spectrum is a mess but if you basically find the average frequency then that will increase with blood flow. I'm not sure you see much pulsing in the capilarys of say the finger where they use this sort of technique.

Graham

p.s. I spent all night thinking about adding vectors!! http://forums.parallax.com/images/smilies/roll.gif

Post Edited (Graham Stabler) : 10/12/2006 9:23:01 AM GMT

Graham Stabler
10-12-2006, 04:28 PM
Just realized I uploaded the old version of the spreadsheet the second time, I'll correct that when I get home.

Graham

daniel
10-12-2006, 06:42 PM
Tracy Allen said...
Chip, regarding the polar addition of vectors. There might be a very constrained way to go about it. Here is a diagram that shows two vectors at angles A1 and A2, lengths L1 and L2,.
I had started doodling about the next step--that of relaxing the equal-length constraint.· Can the vectors be scaled to a unit-length, and their angles "adjusted" by a table lookup (or trig·calculations)·to find the intersection of the original vector thru the circle of unit-length radius?· Then just sum the adjusted angles.

But then it started to look like "slumming" again.

Daniel

daniel
10-12-2006, 07:06 PM
Err, I said that wrong from what I was thinking last night. I'm taking the children on a field trip today, but I'll "re-doodle" and see if it makes sense again.

Daniel

daniel
10-12-2006, 09:31 PM
Tracy Allen said...
If there is a constraint on the lengths, if for example they are both the same, or they are unit vectors, or they have some fixed ratio, then it is possible to solve it geometrically and the lengths drop out.
Where my thinking was going·was·that the GCD of the two lengths seem to me to be a scaling factor, leaving the other factors as the fixed ratio.· If one could quickly find the other factor, then applying Tracy's geometric solution for the angle, and Paul's solution for the resultant length.· Apply this iteratively to each additional polar vector.

I might be wrong.

Daniel

Graham Stabler
10-12-2006, 10:08 PM
I have now pretty much got my algorithm working, I think it had been working for a while actually without me realizing, one problem was my choice of test angle. 90 degrees seemed like a "simple" choice but no, as Xi and Yi are represented as 32 bit fractions they couldn't go negative and so trying to find the solution around 0 or 90 degrees created errors. So I now have an algorithm that will work between 0 and 90 degrees but not at 0 or 90 degrees.

At this point I'm not too sure what to do in terms of the algorithm, as long as I can produce sin/cos in this region then I can produce it for 0-360 degrees by simply making the answer negative or swapping Xi and Yi etc as if you have cos and sin between 0-90 you have everything you need. However if I was to try angles close to 0 or 90 it might also mess up.

As I want to find the x and y components of a vector of length R I'm wondering if ditching the fractional Xi and Yi and going for R*K for the starting vector length might be a better choice. That would allow the variable to be signed.

The other though was to have the algorithm produce the 32bit fraction that represents Sin and Cos and a sign variable, I'd have to keep track of the sign variable myself during the process though.

Thinking about it I think the option of using R*K as the starting length is probably a win-win situation as the only multiplication is a fixed one that has been optimized by chip.

Comments appreciated, I'll be publishing the code so anyone who thinks they might use it should chime in, also wondering how best to present the code, insertable asm "function" or an object.

Graham

Graham Stabler
10-12-2006, 10:40 PM
OK, so it is now coded up to have an input angle and vector length and it works, however I can now understand comments made about the size of numbers. When using the fraction the numbers only got small at the end of the function but the dX and dY soon get small now unless the starting value of R is large. I'm guessing I could take the input R and shift the bits to make it bigger and then shift the resulting value back again.

More playing to do.

Graham

Graham Stabler
10-12-2006, 11:50 PM
It works, I just shift the input R until bit 31 goes high, record the number of shifts needed and then shift back at the end.

The input vector still needs to be of a decent size but in my case I can decide on the size.

I'll post code and example later, I'd still like advice on the best way to do the object.

Graham

Tracy Allen
10-13-2006, 12:24 AM
Peter Verkaik said...
Tracy,
for equal lengths, shouldn't that be A3 = (A1+A2)/2
after all, adding (1,0) to (0,1) yields (1,1)
(0 + pi/2)/2 = pi/4

regards peter


Thanks for catching that, Peter. You must really be awake, while we are burning the midnight oil over here! I like your disproof by counterexample. There is a lot to look at this morning, with the commentary and and Graham's new results. The midnight oil was bubblng.

The reasoning is as follows:

http://forums.parallax.com/attachment.php?attachmentid=43590

for angle:
a4 = A1 chord through parallel lines
a5 = A2 - a4 = A2 - A1 composite angle
a6 = 180 - a5 = 180 - (A2 - A1) = 180 + A1 - A2 composite angle
a7 + a8 = 180 - a6 = 180 - (180 + A1 - A2) = A2 - A1 sum of angles of triangle
2 * a8 = A2 - A1 isocoles triangle, L1=L2, a7=a8
a8 = (A2 - A1)/2 algebra
A3 = A1 + a8 = A1 + (A2 -A1)/2 = (A1 + A2) / 2
therefore,
A3 = (A1 + A2)/2

for length:
a6 = 180 + A1 - A2
construct perpendicular to L3 bisecting a6 and bisecting L3
then L3/2 = L1 * sin(a6/2) = L1 sin(90 + (A1 - A2)/2)
under the condition that a6<=180 degrees.

That should jive with the formula you are discussing with Paul, but I don't see the connection right away. It is a special case.

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

Post Edited (Tracy Allen) : 10/12/2006 6:15:23 PM GMT

cgracey
10-13-2006, 12:48 AM
Graham Stabler said...

p.s. I spent all night thinking about adding vectors!! http://forums.parallax.com/images/smilies/roll.gif

I think we all did! Paul, Peter, you, and Tracy and I have all got the bug now. I had a hard time concentrating on what I was working on because I was tempted to start working on this (as in, "again",·since I never get too far before getting re-overwhelmed). If we can get something figured out, I would love to build polar vector summing hardware into the next generation of the chip. It would make DSP so simple to program and understand. There will at least be a CORDIC system in the hub which all cogs could use. But a polar vector summer would be the cat's meow (or the bee's knee's).

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


Chip Gracey
Parallax, Inc.

Peter Verkaik
10-13-2006, 01:57 AM
Chip, Tracy,

I made a diagram to visualize the iteration process I mentioned.
As it turns out, the vectors that need to be added, move towards
each other until they both have the same length and angle.
At that point the resulting length is just the double value, the angle is the resulting
angle.

In the diagram, start with P0 and Q0. I have drawn circles to find the vector
with smallest length to add to smaller vector. The process is then repeated
with the resulting vector and the remaining vector of the larger vector.
I hope·you understand what I mean.
The iteration should stop when both vectors have·the same angle.

regards peter
·

Post Edited (Peter Verkaik) : 10/12/2006 7:50:36 PM GMT

cgracey
10-13-2006, 02:16 AM
Peter Verkaik said...
Chip, Tracy,

I made a diagram to visualize the iteration process I mentioned.
As it turns out, the vectors that need to be added, move towards
each other until they both have the same length and angle.
At that point the resulting length is just the double value, the angle is the resulting
angle.

In the diagram, start with P0 and Q0. I have drawn circles to find the vector
with smallest length to add to smaller vector. The process is then repeated
with the resulting vector and the remaining vector of the larger vector.
I hope·you understand what I mean.
The iteration should stop when both angles have the same angle.

regards peter
·

Peter,

This seems like a viable concept that would reduce to shifts and adds. This is exciting! I don't understand all the details of what you're doing, but it definitely feels like the right kind of solution. I think you're onto it!!!

Here's Peter's diagram from his last posting:

http://forums.parallax.com/attachment.php?attachmentid=43594

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


Chip Gracey
Parallax, Inc.

Post Edited (Chip Gracey (Parallax)) : 10/12/2006 7:20:57 PM GMT

Peter Verkaik
10-13-2006, 02:27 AM
Chip,
When adding 2 vectors of equal length
the result length = r*sqrt(2 + 2cos(Θa-Θb))
so you would need to find a cos for each iteration.
I can't tell if that is faster than
convert polar to coordinate
add vectors
convert coordinate to polar
as I don't have a clue how fast the angle difference reaches near zero.
Also note the diagram is only for 2 vectors in the same quadrant.

regards peter
·

Tracy Allen
10-13-2006, 02:37 AM
Hi Peter, thanks for the diagram. I'll have to stare at it later, (cuz now I have to focus on my day job). Now it is you burning the midnight oil! You're in Netherlands, right?

A few months ago I had an email exchange about the Given's transform with a person who had seen my write-up and diagrams, and he was investigating the origins of the transforms in the work of the early astronomers. That goes back to the 17th century and all the way back to Ptolemy. Sometimes those early methods have the seed of something new. This person referred me to a diagram he had published on the 3 dimensional transformation equations, noting that the 2 dimensional equations drop right out. Here is that url if anyone is interested...
www.thediagram.com/4_6/moyer.html (http://www.thediagram.com/4_6/moyer.html)

Daniel, I'll be interested to see what you come up with. If the ratio of lengths is constant, it should come in as a multiplier. The CORDIC machinery can do straight multiplcaiton and division and square root, in addition to the fancier functions. They all compute in the same amount of time, so I don't see any advantage in getting rid of a trig function, for example, if it is only to be traded for a multiplcation. _Unless_ the multiplier is built into the hardware! And Chip, that would be fantastic to have the CORDIC engine built into the core. Even a "tranditional", CORDIC engine, very professional, awesome.

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

Peter Verkaik
10-13-2006, 03:09 AM
Hi Tracy,
Yes I live in the Netherlands and it is now 22.10 local time, so not quite midnight here.

Let me try to explain the diagram a little better.
We start with P0 and Q0.
We want to add only vectors of equal length.
Because |Q0| < |P0| I used a circle to find the interception (P1) on vector P0
that has the same length as vector Q0.
R1 is the sum vector Q0+P1
The vector P1 must be subtracted from vector P0. I did this by using an identical
circle (with center P0) to find the interception (P2) that is the remaining vector.
So now we have 2 new vectors: R1 and P2 (R0 = Q0 + P0 = R1 + P2)
Repeat the process for vectors R1 and P2.

Initially I thought the remaining vector length would decrease to zero, but the
diagram shows it is the angle difference between the two new vectors that
decreases to zero.

regards peter

Phil Pilgrim (PhiPi)
10-13-2006, 01:07 PM
Update: Here's a corrected version of this posting, with the erroneous equation for C (i.e. C = (A + B) / 2) removed and replaced by correct formulae, along with a somewhat different conclusion. Rather than leaving the original posting cluttered with overstrikes and edits, I simply fixed it. It's not as proper doing it that way, but it's a lot cleaner!

I'd like to work through a possible alternative to pure polar coordinates to see if it might confer any advantages. In a way, this is a cross between polar and Cartesian coordinates. In this scheme, instead of keeping track of the angle and radius for a vector vA (call them A and rA, respectively), we'll keep track of the sine and cosine of A instead (call them sA and cA, respectively), along with rA, i.e.:

····vA = (rA, sA, cA) = (rA, sin A, cos A)

For signal processing, one will likely need the sine and cosine of the result anyway, so we might as well bring 'em along for the ride from the get-go and, with any luck, avoid having to compute any further transcendental functions.

Consider the following diagram:

····http://forums.parallax.com/attachment.php?attachmentid=43596

We recognize the following identities from the diagram and from trigonometry:

····rC sC = rA sA + rB sB (i.e. yC = yA + yB)
····rC cC = rA cA + rB cB (i.e. xC = xA + xB)
····rC2 = rA2 + rB2 - 2 rA rB cos D (Law of cosines.)
····D = pi - (B - A) = pi + (A - B)

····sin(x + y) = sin x cos y + cos x sin y
····cos(x + y) = cos x cos y - sin x sin y

From these we can derive:

····cos D = cos[pi + (A - B)]
··············= cos pi cos(A - B) - sin pi sin (A - B)
··············= - cos(A - B)
··············= - cos A cos(-B) + sin A sin(-B)
··············= - cos A cos B - sin A sin B
··············= - cA cB - sA sB

So

····rC2 = rA2 + rB2 + 2 rA rB (cA cB + sA sB), leaving

····rC = [rA2 + rB2 + 2 rA rB (cA cB + sA sB)]½, from which we can then compute

····sC = (rA sA + rB sB) / rC
····cC = (rA cA + rB cB) / rC

Now, to summarize what we've got: Given two vectors vA = (rA, sA, cA), and vB = (rB, sB, cB), the sum vC = vA + vB can be computed from the above formulae for rC, sC, and cC. These formulae involve no transcendental functions. They do involve a square root, however, which can be computed easily and converges quickly using Newton's method.

As an addendum, for cycling through the phases of any vector by adding an incremental angle d, first compute sd = sin d and cd = cos d, then use the identities for angle addition:

····r(A + d) = rA
····s(A + d) = sA cd + cA sd
····c(A + d) = cA cd - sA sd

The big question after all this fol de rol is whether it's any better than bopping in and out of Cartesian coordinates. I'm afraid the answer is no. Consider the following, much simpler formulae:

····rC = (xC2 + yC2)½
·········= ((xA + xB)2 + (yA + yB)2)½

····sC = yC / rC
·········= (yA + yB) / rC

····cC = xC / rC
·········= (xA + xB) / rC

Clearly, keeping the x's and y's around confers a significant computational advantage. And, truth be known, in the "modified polar" system derived above, we never really escaped the Cartesian domain, but simply disguised it in trigonometric clothing. Anyway, if was fun working through it, and I hope I didn't waste anyone's time. http://forums.parallax.com/images/smilies/smile.gif

-Phil

Post Edited (Phil Pilgrim (PhiPi)) : 10/14/2006 7:56:28 PM GMT

Graham Stabler
10-13-2006, 04:56 PM
Chip, I wanted to ask something. If an "addpolar" function was created that took two vectors and produced the resultant in polar form would it really matter what was happening behind the scenes?

Say you have a vector, this was created by a vector sum and that sum was done the traditional way by decompositon using cordics, being clever you kept hold of the resultant x,y co-ordinates too. You now want to add a vector, that is just two cordics, one to find the components of the new vector and a second to find the angle and length of the resultant. The new co-ordinates of the resultant is also known.

Those two cordics could be done sequentially or put in the same loop, in the same loop would save a little time because the instructions for t1 and the loop lookup would be common (saves 120 instructions).

So the idea is to have best of both worlds in a way, fairly quick additions and instant rotations and scales. Of course there is the problem that if you have scaled it the x,y's are no longer correct so in that instance you would need to do three cordics to do an addition.

I'm struggling to work out if this is going to be of any use, if you do all your additions and scales seperately then it seems pretty good but if you process interleaves scales and additions then it might be poor, then again the CORDIC is pretty fast so three cordics per rotation/additon pair might not be so bad.

Graham

Peter Verkaik
10-13-2006, 05:19 PM
Phil,
You assume C = (A+B)/2
but this is only true if vectors A and B have the same length.
Adding (1,0) and (0,2) yields (1,2) which does not have a pi/4 angle.
C = (A+B)/2 calculates to pi/4 for vectors along the x and y axis, regardless what vector lengths.
For arbitrary lengths, the calculation of C is more complex.

regards peter

Peter Verkaik
10-13-2006, 06:25 PM
Graham,
If all vector parameters are kept, x, y, r and A,
then when adding 2 vectors you could use x and y to find resulting x and y, then do 1 cordic to find
resulting r and A.
When rotating a vector, simply adjust A, and do 1 cordic to find the new x and y.
When scaling a vector, simply multiply x,y and r or only multiply r and use 1 cordic to find new x and y.
This may well be the fastest solution.

regards peter

Graham Stabler
10-13-2006, 07:56 PM
Peter Verkaik said...
Graham,
If all vector parameters are kept, x, y, r and A,
then when adding 2 vectors you could use x and y to find resulting x and y, then do 1 cordic to find
resulting r and A.
When rotating a vector, simply adjust A, and do 1 cordic to find the new x and y.
When scaling a vector, simply multiply x,y and r or only multiply r and use 1 cordic to find new x and y.
This may well be the fastest solution.
regards peter


I was assumining the vector to be added was in polar co-ordinates but there is no reason why that should not be in x,y,r,A already as you say.

Graham

Graham Stabler
10-13-2006, 07:58 PM
Almost forgot, here is the latest and correctest version of the CORDIC simulator in excel

Phil Pilgrim (PhiPi)
10-13-2006, 10:22 PM
Peter,

Oh cripes! That same realization got me out of bed this morning. What was I thinking?! I had blindly copied the formula from one of Tracy's postings without even examining its premises. In a way, I'm glad. That half-angle formula annoyed the heck out of me. I'll have to rework it, but I'm not at all optimistic that anything can be gained by it. I still believe, however that maintaining the sine and cosine of the angle as part of each vector is more useful than carrying around the angle itself.

-Phil

Paul Baker
10-14-2006, 12:02 AM
Phil, I was thinking a few days ago along the same lines when I was examining the representation of an angle via an exponent to use the power as the notation used, but I didn't delve into the idea very deeply.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Paul Baker (mailto:pbaker@parallax.com)
Propeller Applications Engineer
[/url][url=http://www.parallax.com] (http://www.parallax.com)
Parallax, Inc. (http://www.parallax.com)

Peter Verkaik
10-14-2006, 12:14 AM
Now here is some interesting reading
http://dspace.lib.utexas.edu/bitstream/2152/32/1/arbaughj20424.pdf
Here table lookups are performed to find sin and cos for the first 1/3 of the iterations
What's better, it includes c++ sources to calculate those tables.
Check out chapters 4, 5 and 6.
regards peter

cgracey
10-14-2006, 01:14 AM
Graham Stabler said...
Chip, I wanted to ask something. If an "addpolar" function was created that took two vectors and produced the resultant in polar form would it really matter what was happening behind the scenes?

Say you have a vector, this was created by a vector sum and that sum was done the traditional way by decompositon using cordics, being clever you kept hold of the resultant x,y co-ordinates too. You now want to add a vector, that is just two cordics, one to find the components of the new vector and a second to find the angle and length of the resultant. The new co-ordinates of the resultant is also known.

Those two cordics could be done sequentially or put in the same loop, in the same loop would save a little time because the instructions for t1 and the loop lookup would be common (saves 120 instructions).

So the idea is to have best of both worlds in a way, fairly quick additions and instant rotations and scales. Of course there is the problem that if you have scaled it the x,y's are no longer correct so in that instance you would need to do three cordics to do an addition.

I'm struggling to work out if this is going to be of any use, if you do all your additions and scales seperately then it seems pretty good but if you process interleaves scales and additions then it might be poor, then again the CORDIC is pretty fast so three cordics per rotation/additon pair might not be so bad.

Graham
Graham,

It may, indeed,·turn out that the most practical thing to do is maintain vectors in terms of (X,Y) and use CORDIC to rotate them (with built-in prescaling, to keep·usage as simple as possible). This way, vector adds are still straightforward (two ADD instructions - one for the X's and one for the Y's) and rotation is simple (either a CORDIC or hub instruction). If there was some computationally advantageous technique for adding polar vectors directly, it would be superior since a rotation would amount to a single ADD instruction (for the angle component) and the length would always be immediately accessible. Also, it would allow the programmer to think more clearly about what he's programming, as the polar domain is simpler to contemplate and allows insights and intuitive leaps which are obscured in the cartesian domain.

In a sense, yes, it doesn't matter what transpires in the background, as long as the result is achieved in a timely way, but some novel approach could greatly reduce the circuit complexity which also reduces power and silicon size/cost. That's what intrigues me. I have to believe there is some way to achieve this, even though none of the DSP vendors seem to have come up with anything. It's quite possible that they haven't bothered, since the industry seems quite sated by glorified MACs.

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


Chip Gracey
Parallax, Inc.

Paul Baker
10-14-2006, 03:37 PM
Chip, I found the function to draw a straight horizontal line in polar space as we talked about at dinner w/ Tracy. Its the cosecant (csc) function (and the secant draws a vertical line in polar space). The distance from the origin at Θ=90 degrees is the scaling factor of the function, so 2csc(x) is a horizontal line that is 2 away from the origin at 90 degrees. I think this can be used in your coordiante rotation vector summing concept. I haven't quite worked out how the formula works exactly, but it's there and hopefully will be little more than a single trig function (maybe will require using the secant too)·and an addition. I'll look at this idea some more over the weekend.

Even if its two trig functions and an add, it should be more computationally efficient than converting to cartesian and back.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Paul Baker (mailto:pbaker@parallax.com)
Propeller Applications Engineer
[/url][url=http://www.parallax.com] (http://www.parallax.com)
Parallax, Inc. (http://www.parallax.com)

Graham Stabler
10-14-2006, 06:24 PM
I'm a bit scared to post this because I have probably overlooked something pretty major and I've only considered one case but I think the attached shows a way of adding two polar vectors with only two cordics and that beats normal polar-cartesian-polar by one cordic. If a cordic engine was built into the propeller then you would be cooking on gas. If the method is correct then I am afraid it is a lot cleverer than I am.

Graham

Click other link below for full sized version

http://forums.parallax.com/attachment.php?attachmentid=43627

Post Edited (Graham Stabler) : 10/14/2006 11:28:09 AM GMT

Graham Stabler
10-14-2006, 06:34 PM
How it works in words

The first step is to use my new cordic to find the components of vector2 (this finds both on one go) in the direction of and perpendicular to vector1. This produces a carteasean vector representing the resultant but rotated by -theta1. The second cordic converts this carteasean vector back to polar producing an angle and length. This just needs to be rotated by theta1 which in the polar space is just an addition of theta.

In a sense this method avoids one of the CORDICs that would have converted R1 to carteasean by making the conversion trivial, by simply putting vector1 on the "new" x-axis.

Graham

Peter Verkaik
10-14-2006, 09:42 PM
Graham,
if both vectors also carried x and y, you can leave out the first cordic.
You trade memory for speed. If speed is what you're after, that's
how I would do it. The extra memory required (two 32bit words for x and y)
is worthwile I think to almost halve the time required for the addition of 2 vectors.

regards peter

Graham Stabler
10-14-2006, 11:26 PM
I think they would be carrying the wrong x and y's. You have to find the x,y of vector two with resepect to vector one because the angle of the second vector is with respect to the x-axis.

Graham

Peter Verkaik
10-14-2006, 11:43 PM
Graham,
What I mean is that EVERY vector should have properties x,y,r,A and that these are valid
upon entry and exit of a vector (library) function.

The relation between the properties is
r = sqrt(x*x + y*y)
A = arctan(y/x)

Having x,y,r,A instantly available also allows·tan,sin,cos to be quickly calculated
tan(A) = y/x
sin(A) = y/r
cos(A) = x/r

The add function would then do:
··add the x's and y's
··do 1 cordic to calculate·r and A for the new x and y

The rotate function would do:
· add angle to A
· do 1 cordic to calculate new x and y

The scale function would do:
· multiply r by scale factor
· do 1 cordic to find new x and y


regards peter

Graham Stabler
10-14-2006, 11:53 PM
Sorry I thought your comment was related directly to my previous post rather than the previous discussion.

More likely your polar sums would be:

Find the missing co-ordinates of the vector to be added
Do the add or whatever
calculate the new x,y or r and A

Still two cordics.

It would depend on what you were doing of course, if you were always adding a fixed set of vectors that would be different or in some cases you might flip between systems (for long streams of additions or rotations for example).

Graham

Peter Verkaik
10-14-2006, 11:57 PM
Graham,
You say
Find the missing co-ordinates of the vector to be added
Do the add or whatever
calculate the new x,y or r and A

The first line is not true, because the vector to add already has valid x and y in its properties.
So 1 cordic is all that is required.

regards peter

kelvin james
10-15-2006, 12:26 AM
Not that i really want to get involved here as this is 5 miles over my head. I had previously read about Bresenhams' algorithim for plotting without using trig and uses integer only math. I don't know if this would be applicable to this situation or not, or if it could be manipulated into vector addition. Probably way off base, but thought it might be worth a look. Here is a explanation how it works.

kelvin

www.mandelbrot-dazibao.com/Bresen/Bresen.htm (http://www.mandelbrot-dazibao.com/Bresen/Bresen.htm )

Graham Stabler
10-15-2006, 12:55 AM
Peter, where did those properties come from if not from a cordic? A vector has to start somewhere.

You choose from doing a 1 cordic on every operation or you do what I suggest and use 2 cordics for the additions but only the additions, my method cuts the time for a polar addition down by a third and if as I suspect most DSP consists of alternating scales/rotates and additions then it is no slower than storing all variables.

Graham

p.s. Very interesting Kelvin whether applicable or not!

Tracy Allen
10-15-2006, 01:14 AM
I like where this is going with storing the intermediate results for access on the next cycle. A first step into pipelining, I guess.

Its always fun to visit Parallax, draw diagrams mentally and on napkins. And hear (literally) the progress Chip has made on the speech synth that he posted a few weeks ago. It is clear why he wants the fastest possbible CORDIC routines.

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

Peter Verkaik
10-15-2006, 01:29 AM
Graham,
All properties are set when you define/declare a vector, just as any variable.
Normally you only have properties x and y (or only r and A)
Having only x and y requires no cordic for adding but 2 cordics for rotating.
Having only r and A requires no cordic for rotation but 3 cordics for addition (or 2 in your approach).
Having all 4 properties requires 1 cordic for addition and 1 cordic for rotation.

To illustrate, using some C code
in C a vector would be a structure

struct {
· float x;
· float y;
· float r;
· float A;
} vector;

void setCartesian(vector *v, float x, float y) {
· *v.x = x;
· *v.y = y;
· convertToPolar(v); //make sure r and A are valid
}

void setPolar(vector *v, float r, float A) {
· *v.r = r;
· *v.A = A;
· convertToCartesian(v); //make sure x and y are valid
}

void add(vector *r, vector *a, vector *b) {
· *r.x·= *a.x + *b.x;
· *r.y = *a.y + *b.y;
· convertToPolar(r);
}

void rotate(vector *r, vector *v, float angle) {
· *r.A = *v.A·+ angle;
· *r.r = *v.r;
· convertToCartesian(r);
}

void convertToPolar(vector *v) {
· //cordic code to calculate r and A·using x and y
}

void convertToCartesian(vector *v) {
· //cordic code to calculate x and y·using r and A
}

regards peter

Graham Stabler
10-15-2006, 01:50 AM
It is likely that the program rather than the programmer will define the 4 variables of the vector, this will probably entale a formula or process to get one set of variables and a cordic to get the other, then you do a rotate, scale or addition and you need to do another cordic, that's two cordics for ANY two vector process and potentially 3 depending on how the vectors are created (defined or calculated).

It would very much depend on the processing you intended to do I suspect.

Cheers,

Graham

cgracey
10-15-2006, 01:59 AM
Peter Verkaik said...


Having only x and y requires no cordic for adding but 2 cordics for rotating.

Doesn't this just require 1 cordic? Both X and Y get rotated in a single cordic.

Having all 4 properties requires 1 cordic for addition and 1 cordic for rotation.


I think this is a really interesting idea - keep both forms around all the time. I worry that some error between the representations would build up over time, though. What do you think?

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


Chip Gracey
Parallax, Inc.

Peter Verkaik
10-15-2006, 02:13 AM
Graham,
"It is likely that the program rather than the programmer will define the 4 variables of the vector"
This doesn't matter as long as the library functions are used.

"that's two cordics for ANY two vector process "
That's true, but when having 4 properties, one is in the set function to initialize.
A running vector that is used throughout the program, is usually only initialized once.

Chip,
Having only x and y, to add no cordic is required since there is no need, at that time,
for the length and angle.
1 cordic for addition and 1 cordic for rotation also makes the time required
for the functions more constant I think, which may proof useful in some applications.
If you were to do only additions, and no interest in length and angle, then I would have only x and y.

regards peter

Phillip Y.
10-15-2006, 03:18 AM
I reviewed the link for Bresenhams'·algorithm for plotting.
"These algorithms are called INCREMENTAL, because the position of the next pixel is calculated on the basis of the last plotted one".

In the CNC world this is called interpolation it can ber linear, circular, log, exponential etc, the output·is + and - pulses for·X and·Y servos.
Early Yaskawa/Ysnac CNC machines used an "Interpolation Processor" for this, in 1983·the KM37601·was sold separately by TOKO.
I have·a chip and 75 page user manual with detailed examples of the interpolation calculations/process.

The interpolator uses some control bits to select the type of interpolation, direction and some 24 bit constants for the line properties,
then it would calculate the next move based on the best tracking error·of all the next point's 8·possibilities, it include 45 degree moves.
The calculations were only add and subtract and·the control bits +1,0,-1 [linear 10 clocks·, other 20 clocks for each interpolation cycle].

Another closely related thing is a Velocity Profile Generator to perform acceleration and deceleration of the interpolated steps.
Profiles could be linear, exponential and perhaps Scurve for smoothest starting and stopping.
This is really just another form of interpolation·where one axis is motion and the other axis·is speed or time between steps.

I recommend that·a·FUNCTION ENGINE be made general enough to do CORDIC, Interpolation and many other functions.

Phil Pilgrim (PhiPi)
10-15-2006, 03:27 AM
In a prior posting (http://forums.parallax.com/showthread.php?p=609709), which I've now corrected, I explored the possibility of keeping the length along with the sine and cosine of the angle (instead of the angle itself) of each vector with which to do all the computations. It worked without requiring any further computation of transcendental functions. But I concluded there that it conferred no advantages over a purely Cartesian approach and, in fact was little more than a Cartesian formulation in trigonometric trappings.

Also, I've recently encountered some papers about "rational trigonometry", which you can read about here: en.wikipedia.org/wiki/Rational_trigonometry (http://en.wikipedia.org/wiki/Rational_trigonometry). Mentioning it in this thread may be a red herring, since I doubt its applicability in signal processing apps. But who knows? Perhaps, with a little massaging, it might prove useful. Anyway, it's being touted for its simplification of trig by replacing circular functions and their associated computational difficulties with metrics called "quadrance" and "spread". It's the spread function that I find troubling for signal-processing apps, since the quadrant information gets completely buried (which, the author claims, is okay for plane geometry). Anyway, 'just thought I'd toss it into the fray...

-Phil

Peter Verkaik
10-15-2006, 04:43 AM
Phil,
Thanks for updating your other post. It really shows you can't really do without x and y.
Having also the r and A as vector properties makes both addition and rotation require
only 1 cordic.

The document
http://dspace.lib.utexas.edu/bitstream/2152/32/1/arbaughj20424.pdf
chapter 5 discusses lookup tables to perform the first third of cordic iterations
in a single step. If I understood it correctly, multiple tables can be used
to bring the number of iterations down more.

Chip,
there is also a discussion (chapter 5 and·6) in that document on how to implement that on silicon,
where one needs a fast way to index the lookuptables.
If I find the time I might try·this out·on the javelin, as I already
have a cordic class for it. Just to see how much time can be saved by using tables.

regards peter

Graham Stabler
10-15-2006, 07:03 AM
Peter Verkaik said...


"that's two cordics for ANY two vector process "
That's true, but when having 4 properties, one is in the set function to initialize.
A running vector that is used throughout the program, is usually only initialized once.




Imagine this

program starts
initialize vector1
initialize vector2
repeat
{
calculate vector2 (one cordic to get 4 variables)
add vector2 to vector1 (one cordic to get 4 variables)
}

In this situation you still end up with two cordics, one to redefine the second vector and one to redefine the new resultant.


Peter Verkaik said...

Chip,
Having only x and y, to add no cordic is required since there is no need, at that time,
for the length and angle.



I think chip was just pointing out that it only takes one cordic to rotate an x,y pair (you don't have to convert to polar and back)

I had a bit of a look at the pdf, it looks as though the data tables might need to be quite large?

Graham

Peter Verkaik
10-15-2006, 07:42 AM
Graham,
using the code example functions I described

program starts
setPolar(...); //initialize vector1
setPolar(...); //initialize vector2
repeat
{
add(...); //add vector2 to vector1 (one cordic to get 4 variables)
}

There is no need to calculate vector2 again inside the repeat
as the functions make sure ALL properties are valid all the time,
so you can use them directly when needed.

About the tables,
as I understand it multiple tables can be used that
significantly reduces the total amount of storage required.
According to figure 5.6 (page 65) when using 4 tables
and 32bits coefficients, only about 10 kilobyte (rom) space
is required.

regards peter

Graham Stabler
10-15-2006, 07:56 AM
Peter, your example doesn't do anything useful it justs adds vector2 over and over. Vector 2 might need to be modified on every loop, that might mean completely reloading the x,y values (for example in a position tracking system the values might come from an accelerometer) and that means another cordic to get the new r,A. Or vice versa.

I don't know how to make it any clearer.

Graham

Peter Verkaik
10-15-2006, 08:38 AM
Graham,
The original question was to add 2 vectors in polar properties
and return the result in polar properties.
It turns out we still require x and y to do that.
If you keep vectors in polar·format only you would
· convert vector1 to cartesian
· convert vector2 to cartesian
· add vectors
· convert result to polar
If you keep all properties you would
· add vectors

My point is, the adding itself requires less cordics when having all properties.
The initialization of a vector however also requires a cordic,
so in the end, the total number of cordics in the program
is probably only slighty less. But if you need the components (x and y)
later on, they are instantly available, as are the length and angle.

regards peter

cgracey
10-15-2006, 08:53 AM
Peter Verkaik said...
Chip,
there is also a discussion (chapter 5 and·6) in that document on how to implement that on silicon,
where one needs a fast way to index the lookuptables.
If I find the time I might try·this out·on the javelin, as I already
have a cordic class for it. Just to see how much time can be saved by using tables.

regards peter


Peter (and Graham),

I got up to page 70 or so of that document and when I saw that he was proposing 2048K table entries, I quit reading, as this would not be practical in our chip.

I really appreciate all the thought you, Graham, Tracy, Paul·and the others have put into this question. You guys are much better versed in the math than I am, so I have to read what you've written very carefully to absorb it.·Paul Baker and I·had dinner with Tracy Allen last night and we were discussing something like Graham has proposed: rotating a polar pair to sum against the X axis and maybe shortening the work. I must read more of what Graham said regarding this.

The next chip will have a single-cycle signed/unsigned 16x16 multiplier built in, so I'm kind of leaning towards a simple 16-bit (X, Y, and angle) CORDIC system in the hub which would be pipelined. It could rotate (X,Y) by an angle and also rotate Y to 0 to return the length and angle. This would provide all the basic functions that we need and the pre-/post-scaling could be done externally in software using two multiply instructions. This way,·with 14 instructions in-between each CORDIC exchange, you could do all the pre- and post-work. I'm still anxious to know about any possible labor-saving devices, though. I have a hunch that there is some technique for doing this that is way out of the box and only afterwards could it be tied back to the math identities. It's just very hard to get there from here.

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


Chip Gracey
Parallax, Inc.

Tracy Allen
10-15-2006, 11:42 AM
The thing in Jason Arbaugh's paper that really caught my attention was the Hybrid CORDIC algorithm described at the end of chapter 4. Thanks for that reference, Peter. The crux of the method is that the value of tangent(theta) approaches close to theta when theta is small., or, to say the the same thing, theta by itself becomes a better and better approximation to arctan(theta) as theta becomes small. That is what happens in CORDIC successive approximation...
arctan(1 / 2^i) => 1 / 2^i <--- where by "=>" it means within an error that decreases as i increases.

The hybrid algorithm uses the standard arctan radix up through about 1/3 of the iterations, and then (as shown in a paper by Wang et al., Arbaugh's reference 52), the algorithm stops reading from the arctan table and instead uses the value 1 / 2^i directly. That saves having to read from the table for 2/3 of the iterations, but there is more. At that point, with the residual angle now in a staight binary radix, the 1's and 0's are one and the same as the sequence of + and - rotations, so there is no longer any need for test and subtract. The algorithm just has to step through and add 1/2^i when a bit is one, and subtract that when a bit is zero. The last 2/3 could all be done in a special purpose combinatorical adder without any further looping, so it can be "parallel".

The CORDIC gain is not affected, because it has almost converged at the 1/3 point, and the angles continue to accumulate. The reference Albaugh has for the hybrid algorithm is his number 52, (Wang, Piuri and Swartzlander, "Hybrid CORDIC Algorithms" IEEE transactions on Computers, v46, pp1202-1207, 1967).

Arbaugh goes on to develop the table method to decrease the execution time for the first 1/3 of the algorithm, which is then combined with the hybrid method. I don't really have a handle on the table method, but there are significant memory requirements.

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

Post Edited (Tracy Allen) : 10/15/2006 5:58:03 AM GMT

Tracy Allen
10-15-2006, 12:02 PM
Just to make the point stronger about what is going on in the hybrid algorithm, here is a spreadsheet that shows how rapdily arctangent of 1 / (2^i) converges to 1 / (2^i). There is the difference, and also the difference normalized to the arctangent value (the percent diffeerence if you will). This convergence allows the substitution, within some small error bound.




I 1 / 2^I atn(1 / 2^I) difference normalized
0 1.000000000000E+00 7.853981633974E-01 2.146018366026E-01 2.732395447352E-01
1 5.000000000000E-01 4.636476090008E-01 3.635239099919E-02 7.840521614581E-02
2 2.500000000000E-01 2.449786631269E-01 5.021336873136E-03 2.049703761562E-02
3 1.250000000000E-01 1.243549945468E-01 6.450054532386E-04 5.186807780334E-03
4 6.250000000000E-02 6.241880999596E-02 8.119000404265E-05 1.300729764760E-03
5 3.125000000000E-02 3.123983343027E-02 1.016656973172E-05 3.254361056187E-04
6 1.562500000000E-02 1.562372862048E-02 1.271379523169E-06 8.137491082010E-05
7 7.812500000000E-03 7.812341060101E-03 1.589398988889E-07 2.034472095702E-05
8 3.906250000000E-03 3.906230131967E-03 1.986803302824E-08 5.086242324960E-06
9 1.953125000000E-03 1.953122516479E-03 2.483521181242E-09 1.271564461670E-06
10 9.765625000000E-04 9.765621895593E-04 3.104406805406E-10 3.178913579284E-07
11 4.882812500000E-04 4.882812111949E-04 3.880510171007E-11 7.947285461817E-08
12 2.441406250000E-04 2.441406201494E-04 4.850638228755E-12 1.986821457973E-08
13 1.220703125000E-04 1.220703118937E-04 6.063297921469E-13 4.967053681939E-09
14 6.103515625000E-05 6.103515617421E-05 7.579122740650E-14 1.241763471370E-09
15 3.051757812500E-05 3.051757811553E-05 9.473904272845E-15 3.104408953090E-10
16 1.525878906250E-05 1.525878906132E-05 1.184238457622E-15 7.761025156475E-11
17 7.629394531250E-06 7.629394531102E-06 1.480300189610E-16 1.940259064563E-11
18 3.814697265625E-06 3.814697265607E-06 1.850385824924E-17 4.850675416913E-12
19 1.907348632813E-06 1.907348632810E-06 2.313035220715E-18 1.212696609800E-12
20 9.536743164063E-07 9.536743164060E-07 2.891558723689E-19 3.032019080252E-13
21 4.768371582031E-07 4.768371582031E-07 3.615771893592E-20 7.582823258190E-14
22 2.384185791016E-07 2.384185791016E-07 4.526332311890E-21 1.898481372109E-14
23 1.192092895508E-07 1.192092895508E-07 5.691002614365E-22 4.773959005888E-15
24 5.960464477539E-08 5.960464477539E-08 7.279189390467E-23 1.221245327088E-15
25 2.980232238770E-08 2.980232238770E-08 0.000000000000E+00 0.000000000000E+00
26 1.490116119385E-08 1.490116119385E-08 0.000000000000E+00 0.000000000000E+00
27 7.450580596924E-09 7.450580596924E-09 0.000000000000E+00 0.000000000000E+00
28 3.725290298462E-09 3.725290298462E-09 0.000000000000E+00 0.000000000000E+00
29 1.862645149231E-09 1.862645149231E-09 0.000000000000E+00 0.000000000000E+00
30 9.313225746155E-10 9.313225746155E-10 0.000000000000E+00 0.000000000000E+00
31 4.656612873077E-10 4.656612873077E-10 0.000000000000E+00 0.000000000000E+00


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

Post Edited (Tracy Allen) : 10/15/2006 5:54:45 AM GMT

cgracey
10-15-2006, 12:48 PM
Tracy Allen said...
The hybrid algorithm uses the standard arctan radix up through about 1/3 of the iterations, and then (as shown in a paper by Wang et al., Arbaugh's reference 52), the algorithm stops reading from the arctan table and instead uses the value 1 / 2^i directly. That saves having to read from the table for 2/3 of the iterations, but there is more. At that point, with the residual angle now in a staight binary radix, the 1's and 0's are one and the same as the sequence of + and - rotations, so there is no longer any need for test and subtract. The algorithm just has to step through and add 1/2^i when a bit is one, and subtract that when a bit is zero. The last 2/3 could all be done in a special purpose combinatorical adder without any further looping, so it can be "parallel".

The CORDIC gain is not affected, because it has almost converged at the 1/3 point, and the angles continue to accumulate. The reference Albaugh has for the hybrid algorithm is his number 56, (Wang, Piuri and Swartzlander, "Hybrid CORDIC Algorithms" IEEE transactions on Computers, v46, pp1202-1207, 1967).


Arbaugh was talking about double-precision IEEE 754 quality results, too, which have 53 bits of (including the leading 1) of mantissa. For a 16-bit system, we wouldn't even get past his 1/3 point where the adding could get paralleled. Just an observation. This is good to know, though for a high-quality double-precision CORDIC implementation that could be in software.

Thanks for explaining all that Tracy, and showing the convergence. I didn't glean nearly as much as you did from that paper.

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


Chip Gracey
Parallax, Inc.

Post Edited (Chip Gracey (Parallax)) : 10/15/2006 5:51:59 AM GMT

Tracy Allen
10-15-2006, 01:30 PM
Chip, while there was mention of the 53 bit mantissa floating point, I don't see any specific reference to that specific bit depth in connection with the logic behind the hybrid algorithm. The formula relating the sufficient number of iterations using the standard arctan radix came from the formula given by Arbaugh as,
n = (N - log2(3)) / 3
where I presume N is the bit depth. The constant log2(3) = 1.58, so that would put n=14/3 = 5, rounded off. Say you play it on the safe side and take the 16 bit algorithm up to iteration 8 and switch at that point to straight 1/ (2^n) instead of the table lookup of arctan(.).

From the spreadsheet,


I 1 / 2^I atn(1 / 2^I) difference normalized
8 3.906250000000E-03 3.906230131967E-03 1.986803302824E-08 5.086242324960E-06



At that point the approximation of column two to column 3 is one part in 100 million absolute and one part in a million normalized to the "correct" arctan value.

That formula comes from an error analysis that is presumably contained in the paper by Wang et al., and i have not read that. But it seems to me that the argument would hold true at any bit depth. The summary implies that the margin of error in the "1/3 rule" is something like one least significant bit, something like that.

Note that theta is always at least a hair greater than arctan(theta), and that is one condition that ensures that there are no "blind spots" that the successive approximation cannot reach.

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

Post Edited (Tracy Allen) : 10/15/2006 6:36:24 AM GMT

cgracey
10-15-2006, 02:10 PM
Tracy,

Sorry I got some concepts twisted around.

It seems that in order to exploit this halving phenomenon that occurs somewhere down the lookup table, you would not use, say, the value $20000000 as a 45-degree starting point, but rather something a little different that would resolve to simple powers of two from the bottom of the table upwards, right? If you look at my table below, you can see that at some point, the values just start halving. The trouble is, they need to follow the pattern: $80, $40, $20, ...$01 at the end, in order to lend themselves to parallelization, not these non-power-of-two values that have more than one bit set in them, causing them to not halve cleanly. This would mean that some angle scaling would always need to be performed, before and after a CORDIC transform, right? Or, some non-power-of-two number would represent 2pi.

···32-bit······ 20-bit

···$20000000··· $20000· ;1
···$12E4051E··· $12E40· ;2
···$09FB385B··· $09FB4· ;3
···$051111D4··· $05111· ;4
···$028B0D43··· $028B1· ;5
···$0145D7E1··· $0145D· ;6
···$00A2F61E··· $00A2F· ;7
···$00517C55··· $00518· ;8
···$0028BE53··· $0028C· ;9
···$00145F2F··· $00146· ;10
···$000A2F98··· $000A3· ;11
···$000517CC··· $00051· ;12
···$00028BE6··· $00029· ;13
···$000145F3··· $00014· ;14
···$0000A2FA··· $0000A· ;15
···$0000517D··· $00005· ;16
···$000028BE··· $00003· ;17
···$0000145F··· $00001· ;18
···$00000A30··· $00001? ;19
···$00000518··· $00000· ;20
···$0000028C··· $00000· ;21
···$00000146··· $00000· ;22
···$000000A3··· $00000· ;23
···$00000051··· $00000· ;24
···$00000029··· $00000· ;25
···$00000014··· $00000· ;26
···$0000000A··· $00000· ;27




Tracy Allen said...
Chip, while there was mention of the 53 bit mantissa floating point, I don't see any specific reference to that specific bit depth in connection with the logic behind the hybrid algorithm. The formula relating the sufficient number of iterations using the standard arctan radix came from the formula given by Arbaugh as,
n = (N - log2(3)) / 3
where I presume N is the bit depth. The constant log2(3) = 1.58, so that would put n=14/3 = 5, rounded off. Say you play it on the safe side and take the 16 bit algorithm up to iteration 8 and switch at that point to straight 1/ (2^n) instead of the table lookup of arctan(.).

From the spreadsheet,


I 1 / 2^I atn(1 / 2^I) difference normalized
8 3.906250000000E-03 3.906230131967E-03 1.986803302824E-08 5.086242324960E-06



At that point the approximation of column two to column 3 is one part in 100 million absolute and one part in a million normalized to the "correct" arctan value.

That formula comes from an error analysis that is presumably contained in the paper by Wang et al., and i have not read that. But it seems to me that the argument would hold true at any bit depth. The summary implies that the margin of error in the "1/3 rule" is something like one least significant bit, something like that.

Note that theta is always at least a hair greater than arctan(theta), and that is one condition that ensures that there are no "blind spots" that the successive approximation cannot reach.


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


Chip Gracey
Parallax, Inc.

Graham Stabler
10-15-2006, 04:58 PM
Tracy, I like the fact the 1/2^i can be used directly but that is only true if your angle is represented in radians, it doesn't really matter I suppose as long as the bit of the table you do use is in radians and you are prepared to accept potential conversions I suppose, ah I think this might be what Chip is saying.

In terms of loosing the decision, on my examples they do seem to converge to adding but only for the last three steps.

Graham

Graham Stabler
10-15-2006, 05:24 PM
Duh, now i get it. We don't actually have 1/2^i to use directly anyway so all you do it use the table of angles (as they are now) and then at some point start using the previous angle from the table divided by two as that equals the next value.

http://forums.parallax.com/attachment.php?attachmentid=43635

Peter Verkaik
10-15-2006, 05:52 PM
Chip,
Pages 58-61 explain how to get easy index for the tables.
The traditional arctangent radix of 2^-i was choosen to
reduce multiplications to simple shifts. (pages 58-59)
The parallel arctangent radix uses an other angle,
where the bits of the angle are grouped and used as index
for the tables. (pages 60-61). This is the tricky part.

I checked my Javelin arctangent table for use with this approach and get
· private final int T = (short)20861; //16384/(pi/4)
· private static int[] atnTable = new int[]{16384,9672,5110,2594,1302,
··········································· 652,326,163,81,41,20,10,5,3,1,0};
································· //T*atn(1/1),T*atn(1/2),T*atn(1/4)...T*atn(1/32768)
· private static int[] atnTable2= new int[]{20861,10430,5215,2608,1304,
············································652,32 6,163,81,41,20,10,5,3,1,0};
································· //T*(1/1),T*(1/2),T*(1/4)...T*(1/32768)

And indeed, from 1/32 on (one third of the table), the values are identical.

regards peter

Graham Stabler
10-15-2006, 06:30 PM
Getting the index is one thing but the ROM would take up a lot of silicon which is what I think put chip off.

Graham

Peter Verkaik
10-15-2006, 07:09 PM
On page 64 is a formulae to calculate storage requirements.
For n=16 bits and t=4 tables 192 bytes would be required
For n=32 bits and t=8 tables 768 bytes would be required
There is a diagram on page 65.

regards peter

Graham Stabler
10-15-2006, 07:18 PM
Is that bytes per table? So would you actually need 8*768 bytes?

Peter Verkaik
10-15-2006, 07:20 PM
No, total amount.
regards peter

Graham Stabler
10-15-2006, 07:31 PM
What was the 10k you refered to earlier?

Graham Stabler
10-15-2006, 07:35 PM
Chip,

I realized that my simplified explanation of the 2-cordic add probably made it sound worse than it is and my last diagram didn't help much. Here is the better explanation, notice the rotate to the x-axis is really just calculating the angle between the two vectors

click on cordic2.jpg for bigger version

http://forums.parallax.com/attachment.php?attachmentid=43637

Peter Verkaik
10-15-2006, 07:41 PM
That was n=32 and t=4.
The diagram shows it to be almost 10000.
Calculation yields 6144 (the diagram has a log scale, so the middle between 1E+02 and 1E+04
is 1E+03 (1000).

regards peter

Graham Stabler
10-15-2006, 07:57 PM
I think we will have to wait for Chip's judgement, as cog ram is currently 2k, I'm guessing that even 500 bytes is considered quite a lot but its certainly better than 2k or 10k.

I just realized something else, if you know x,y,r,A for JUST ONE of the two vectors in the method above you can workout the vector sum with one cordic. It doesn't matter which one you just change the order of the addition. This could be VERY handy for certain things like vector intergrations, you keep all for variables of the main vector and just add a constantly rotating and scaling vector too it. WRONG!! (I thought it gave you the new x,y of the main vector but it doesn't.

Graham

Post Edited (Graham Stabler) : 10/15/2006 1:08:26 PM GMT

Tracy Allen
10-16-2006, 01:14 AM
Graham, I think you are right, that instead of reading values from the table for i=n to 15, you take the n-1th value from the table and successivey divide it by two for the successive rotations, and add or subtract, depending on the bit in the final residual angle being 1 or 0. That is for the case of when the test is done on the accumulated angle to determine sine and cosine of theta. It is different when determining driving y to zero, where the test is done on y, not angle (determining arctangent from x and y).

Alternatively the binary value to represent pi/4 could actually be the binary for pi/4 = %1011010100000101 (16 bits). Then the 1 / 2^i values could be used directly. But there would be a scaling required to get brads and degrees, and I think the math to calculate for different quadrants might get pretty messy.

This algorithm and the Arbaugh table algorithm mesh together. The procedure Peter is talking with the tables covers the first 1/3 of the calculation, and the Wang hybrid approximation covers the final 2/3 of the way to the answer.

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

Peter Verkaik
10-26-2006, 03:54 AM
Another interesting approach.
http://www.ece.utexas.edu/~adnan/comm-06/cordic.pdf

regards peter

Wulffy
04-18-2007, 12:22 PM
Peter Verkaik said...
Another interesting approach.
http://www.ece.utexas.edu/~adnan/comm-06/cordic.pdf

regards peter

http://forums.parallax.com/images/smilies/mad.gif· That link is broken - does anyone have the .pdf saved, that they can legally make available?
Also, one quick question:·Is·the concept of spherical·trigonometry covered herein at all?· I understand that the concepts presented herein lend themselves nicely to planar geometry, however my application (great circle navigation) could surely benefit from the application of spherical trig...
Regardless, I plan on researching this to it fullest, in the interests of getting 'up to speed' as I am embarking on a journey that will either require the use of a uM-FPU, or an understanding/implementation of the concepts presented in this thread (or some darn lucky porting ;-)
In a cursory review of the content presented herein, I feel that I am walking in the land of giants, looking around in sheer awe.· I am humbled by·volume of knowledge that is seemingly offered up with no expectation of reciprocation - kudos and thanks!

Harrison.
04-18-2007, 01:03 PM
The cordic stuff just moved to users.ece.utexas.edu/~adnan/comm-07/cordic.pdf (http://users.ece.utexas.edu/~adnan/comm-07/cordic.pdf) . Professors tend to move things around so that students can't copy stuff from last semester.

Wulffy
04-18-2007, 08:38 PM
Thanks Harrison!

-t

Tracy Allen
04-18-2007, 10:05 PM
Wulffy,

There is a succinct diagram for the rotational transform in spherical coordinateds here:
www.thediagram.com/4_6/moyer.html (http://www.thediagram.com/4_6/moyer.html)
And there is oodles of spherical trig to be found in technical astronomy books. I don't know of any CORDIC shortcuts in relation to these things.

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

Wulffy
04-18-2007, 10:42 PM
Tracy Allen said...
Wulffy,

There is a succinct diagram for the rotational transform in spherical coordinateds here:
www.thediagram.com/4_6/moyer.html (http://www.thediagram.com/4_6/moyer.html)
And there is oodles of spherical trig to be found in technical astronomy books. I don't know of any CORDIC shortcuts in relation to these things.


Mr. Allen,
Thanks so much for taking the time to point out this resource.
In the interests of possibly helping someone else with a similar need, I have also located some pertinent information, as it relates to Great Circle navigation, here:
http://williams.best.vwh.net/avform.htm
Have a wonderful day!
-t

Post Edited (Wulffy) : 4/18/2007 3:52:03 PM GMT

Wulffy
04-20-2007, 10:22 AM
Bean (Hitt Consulting) said...
Here is a document that does a decent job of explaining CORDIC atan/hyp algorithm.

http://www.gmw.com/magnetic_measurements/Sentron/sensors/documents/AN_125KIT_manual.pdf#search=%22an125_kit%20cordic% 22

Bean.
·

Bean - Do you have a copy of this manual archived?· I've dug and dug at GMW's site to no avail.
Please advise.· TIA.
-t

Graham Stabler
04-20-2007, 04:32 PM
Try this:

www.voidware.com/cordic.htm (http://www.voidware.com/cordic.htm)

it has algorithms for most of the cordic routines though it doesn't explain them.

Wulffy
04-21-2007, 12:00 PM
Graham Stabler said...
Try this:

www.voidware.com/cordic.htm (http://www.voidware.com/cordic.htm)

it has algorithms for most of the cordic routines though it doesn't explain them.
Graham,
Thanks for pointing me toward the link.· I am still trying to wrap my mind around the mechanics of the implementations of the·CORDIC algorithm, thus it is actually the theory·that I am after.· However, I so appreciate your input, as I did not make that clear in my initial queries.
Regarding the theories, it is finally starting to sink in ... I think ... I am just looking for a couple more papers to help solidify what·I think is on the tip of my 'mental tongue'.
...
Also, I am trying to port the attached to a 32bit NXP MCU with some simple BASIC on it - integers only, no floats, etc.· It would be·inappropriate for me to request assistance specific to the product, however,·I will ask for assistance with understanding the attached as it relates to·the following two topics:
Struggle #1:
In the attached, I am really struggling with grasping the concept of fixed point notation.· The preamble describes that the number unit is made up of 32 bits.· If I am understanding it correctly, b31 = sign bit, b30 & b29 are the integral portions, and b28 to b0 is the fraction portions.· Now, I fully get the sign and integral portions, but I am struggling with how exactly the use of the remaining 29 bits is able to accurately represent an arbitrary·fractional number, without some obnoxious scaling (with respect to a fixed precision) having to be implemented.· Does the implementation·build the fractional portion·from the inferred decimal point out to the right to b0 (i.e. of the <0·decimal portion, the LSb=b28 and the MSb=b0)?· In re-reading what I just typed, I think I have confused the matter further for myself...·http://forums.parallax.com/images/smilies/rolleyes.gif
Struggle #2:
Also, related to the attached, the 3rd and 4th·lines of code is causing me heartburn (Hold on, wait, don't go ... I think I can explain where you won't think me a total bumbling idiot.?.http://forums.parallax.com/images/smilies/blush.gif ):
Obviously, the code: #define One (010000000000L>>1) is initializing the var One with a single bit set.· What is not obvious to me is how and which bit, exactly, is being set.· The number is 10,000,000,000.· Shifting it right one bit seems to me to serve no logical purpose:
10,000,000,000·equals ·1001010100000010111110010000000000·(>32 bits and not seemingly a good pattern...)
After the shift, it is 0100101010000001011111001000000000, if I am understanding things correctly.
Shouldnt the·number be·· 00100000000000000000000000000000?· If so, what was the logic in using 10,000,000,000?
...
I think that if I can get pointed in the right direction with these two struggles, and get my mind wrapped around the whys and hows, then the 4th line (#define HalfPi (014441766521L>>1))·may make sense to me at that point...
...
As always, thank you, kind sirs!· Your assistance and graciousness to this point is, and will be, greatly appreciated!
-t

Post Edited (Wulffy) : 4/21/2007 5:07:57 AM GMT

Tracy Allen
04-21-2007, 11:36 PM
Wuffy,

Did you also read this thread, Cordic Object (http://forums.parallax.com/showthread.php?p=642142), and also look at the reference there to Turner's book that has code in Matlab.

Often the core algorithm is implemented to cover the domain of convergence only, and then a wrapper is needed that will extend the domain of validity to a range appropriate to a particular application. That way you can get either speed or generality, two qualities that are often in conflict.

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

Graham Stabler
04-22-2007, 03:00 AM
Struggle 1:

Generally in a binary fraction the msb equals 1/2 then the next less significant bit equals 1/4 and so on. One is therefore obtained when all bits are one.

For the whole combination to be one the number should be

00100000_00000000_00000000_00000000

Sign is 0
The integral is one
The fraction is zero.

Struggle 2:

I don't understand what the L represents in the line

#define One (010000000000L>>1)

Nor can I imagine why the shift is used which leads me to suspect all it not as it seems, I suspect the result of this strange operation may be what I suggest in part 1.

To understand the basic concepts I found a mix of wikipedia and Tracy's site good after that having a go even with floats in excel was helpfull. You can see from the link Tracy has just posted the question of "what is one" is key in some cases.

Graham

Wulffy
04-22-2007, 09:20 PM
Tracy Allen said...
Wuffy,

Did you also read this thread, Cordic Object (http://forums.parallax.com/showthread.php?p=642142), and also look at the reference there to Turner's book that has code in Matlab.

Often the core algorithm is implemented to cover the domain of convergence only, and then a wrapper is needed that will extend the domain of validity to a range appropriate to a particular application. That way you can get either speed or generality, two qualities that are often in conflict.


Yeah, I read that thread, and have the chapter referenced therein.· Thanks for pointing that out.· The concept of convergence is escaping me at the moment.· This whole exercise is humbling... :-)

Wulffy
04-22-2007, 10:16 PM
Graham Stabler said...
Struggle 1:

Generally in a binary fraction the msb equals 1/2 then the next less significant bit equals 1/4 and so on. One is therefore obtained when all bits are one.

For the whole combination to be one the number should be

00100000_00000000_00000000_00000000

Sign is 0
The integral is one
The fraction is zero.


I concur with your assessment regarding the sign, integral, and mantissa(fraction).· I am guessing that I will need to continue to be a bit more of a sponge, before I can truly grasp the implementation of the manipulation of the fractional bits.· Regardless, Thanks.

Struggle 2:

I don't understand what the L represents in the line - Neither do I... :) I am suspecting that it is a directive operand that is used to denote a long, but that is speculation.· I have looked at a manual for C as well as the C Preprocessor, to get it set in stone what it means, but it has not yet bubbled to the surface.

#define One (010000000000L>>1)

Nor can I imagine why the shift is used which leads me to suspect all it not as it seems, I suspect the result of this strange operation may be what I suggest in part 1.·- I fully concur regarding your suspicions.· The logic of it escapes me as well, and causes me to be a bit suspect.

To understand the basic concepts I found a mix of wikipedia and Tracy's site good after that having a go even with floats in excel was helpfull. You can see from the link Tracy has just posted the question of "what is one" is key in some cases. - I'll continue to absorb.· I love the concept of having a library that facilitates having a dynamic representation of ONE.· i.e. for those calculations where the·the values being manipulated·are large integers, a 0100 being representative of·ONE allows the majority of the variable's bit space to be representative of the integral, where as a value of 01000000000 would allow for much higher precision (and of course,·smaller integral representations...).

Graham
My approach to date, has been trying to find a library that would lend itself to porting to my platform with minimal understanding required of what is actually transpiring - i.e. the easy way out...
I am finding that·any implementation is actually heavily predicated on a full and finite understanding of the algorithm's mechanics,·having to have a working solution for·handling non-integer·numeric representation in a non-float environment, as well as the limitation on operands that are available in a lower level BASIC implementation.·
Some of the source code out there - i.e. Peter's Java Implementation, or the source I attached, employs concepts specific to environments that I am not familiar with nor operating in.· Not having programmed in an OOP environment before, lends itself to looking at a hunk of code, and initially thinking that it's translation to my environment would be easily accomplished, only to realize that concepts such as recursion, scope, and operator overloading·serve to truly muddy the waters (at least for me... :-).··Not being able to pass/return parameters, but rather having to use·global variables serves to confuse things further as well, during my porting attempts to date.
I think my initial approach·was/is flawed.· I am finally relegating myself to the fact that I must grasp the CORDIC principles first, in detail·(which require an understanding of the underlying trigonometric fundamentals - it has been a lot of years since high-school·:-), and then use the various works available as a frame of reference when drafting my own CORDIC engine...
I·know that when I am am done, I will appreciate what goes into developing same...
At any rate, thank you all for your assistance thus far.· I think I'll give Tracy's code one more shot at porting (yes, I am being lazy...http://forums.parallax.com/images/smilies/rolleyes.gif)
-t

·

Graham Stabler
04-22-2007, 10:59 PM
Well if I can manage to finish my cordic object then you won't need to understand it to use it but an understanding of the number representations on input and output will be required, I was hoping to explain these with many examples in the demo code.

Graham

Peter Verkaik
04-22-2007, 11:42 PM
Graham,
In C, numbers starting with 0, are octal.
A suffix L identifies a long (32bits).
010000000000L>>1 = 001_000_000_000_000_000_000_000_000_000_000>>1 (bin)
················ = 000_100_000_000_000_000_000_000_000_000_000
················ =· 00_100_000_000_000_000_000_000_000_000_000
················ = 0010_0000_0000_0000_0000_0000_0000_0000 (bin)
················ = 2000_0000 (hex)

So you got the number right.

regards peter

Wulffy
04-23-2007, 12:33 AM
Peter Verkaik said...
Graham,
In C, numbers starting with 0, are octal.
...
Doh! http://forums.parallax.com/images/smilies/shocked.gif
Oh my, sometimes we can't see the forest through the trees - this is too funny - Peter, Thank you for showing me the errors of my ways in that... ROTFL http://forums.parallax.com/images/smilies/tongue.gif

Post Edited (Wulffy) : 4/22/2007 5:37:40 PM GMT

Graham Stabler
04-23-2007, 03:53 PM
This demonstrates the problem with search engines, exactly what would I have searched for to find this information "confusing #define with an L and it happens to start with a zero if that makes any difference not that I have any reason to think it does" :)

Thanks Peter!

Graham

Peter Verkaik
04-23-2007, 04:20 PM
Octal numbers in C are one of the lesser attractive syntax.
They (Kernighan and Ritchie) should have used 0o.... for octal, as they use 0x.... for hexadecimal.
Then also 0b.... should be used for binary (some C's actually have this, others use %.... for binary).

The suffix L is a short notation for typecast to long
1L equals (long)1
I guess you just need to know these things, but if you're not into C, it is confusing.

regards peter