Forums

Understanding the CORDIC algorithm with Python

Started by Chico May 6, 2009
Hello, I've been trying to figure out the way the CORDIC algorithm
works and I am not sure if this is the right place to post this, so
feel free to point me in the right direction to get help.

I usually code a quick prototype in python for verification and
understanding of an algorithm I need to implement. For those of you
who don't know Python but know C, you should be able to follow it
since python has a similar syntax to C. I've been struggling with this
CORDIC algorithm, seems straight forward from the wikipedia page, but
somehow I messing something up.


import math

# I know CORDIC is only valid for inputs between
# -pi/2 and pi/2; I am not exactly sure what I need
# to do add to make it where any input is acceptable
# I believe is keep on adding/subtracting pi, but I don't
# get why this works?
def cordic_trig(beta,N=40):
    # in hardware, put this in a table.
    def K_vals(n):
        K = []
        acc = 1.0
        for i in range(0, n):
            acc = acc * (1.0/math.sqrt(1 + 2.0**(-2*i)))
            K.append(acc)
        return K

    #K = K_vals(N)
    K = 0.6072529350088812561694
    # emulation for hardware lookup table
    atans = [math.atan(2.0**(-i)) for i in range(0,N)]
    #print K
    #print atans
    x = 1
    y = 0

    for i in range(0,N):
        d = 1.0
        if beta < 0:
            d = -1.0

        x = (x - (d*(2.0**(-i))*y))
        y = ((d*(2.0**(-i))*x) + y)
        # in hardware put the atan values in a table
        beta = beta - (d*math.atan(2**(-i)))
    return (K*x, K*y)

if __name__ == '__main__':
    beta = math.pi/6.0
    print "Actual cos(%f) = %f" % (beta, math.cos(beta))
    print "Actual sin(%f) = %f" % (beta, math.sin(beta))
    cos_val, sin_val = cordic_trig(beta)
    print "CORDIC cos(%f) = %f" % (beta, cos_val)
    print "CORDIC sin(%f) = %f" % (beta, sin_val)
On Wed, 06 May 2009 14:33:14 -0700, Chico wrote:

> Hello, I've been trying to figure out the way the CORDIC algorithm > works and I am not sure if this is the right place to post this, so > feel free to point me in the right direction to get help. > > I usually code a quick prototype in python for verification and > understanding of an algorithm I need to implement. For those of you > who don't know Python but know C, you should be able to follow it > since python has a similar syntax to C. I've been struggling with this > CORDIC algorithm, seems straight forward from the wikipedia page, but > somehow I messing something up.
> x = (x - (d*(2.0**(-i))*y)) > y = ((d*(2.0**(-i))*x) + y)
You're calculating y[i+1] from x[i+1], when it should be x[i]. Rewriting the above as: (x,y) = (x - (d*(2.0**(-i))*y), (d*(2.0**(-i))*x) + y) to update x and y concurrently gives the correct result.
On May 6, 7:41=A0pm, Nobody <nob...@nowhere.com> wrote:
> On Wed, 06 May 2009 14:33:14 -0700, Chico wrote: > > Hello, I've been trying to figure out the way the CORDIC algorithm > > works and I am not sure if this is the right place to post this, so > > feel free to point me in the right direction to get help. > > > I usually code a quick prototype in python for verification and > > understanding of an algorithm I need to implement. For those of you > > who don't know Python but know C, you should be able to follow it > > since python has a similar syntax to C. I've been struggling with this > > CORDIC algorithm, seems straight forward from the wikipedia page, but > > somehow I messing something up. > > =A0 =A0 =A0 =A0 x =3D (x - (d*(2.0**(-i))*y)) > > =A0 =A0 =A0 =A0 y =3D ((d*(2.0**(-i))*x) + y) > > You're calculating y[i+1] from x[i+1], when it should be x[i]. > > Rewriting the above as: > > =A0 =A0 =A0 =A0 =A0(x,y) =3D (x - (d*(2.0**(-i))*y), (d*(2.0**(-i))*x) + =
y)
> > to update x and y concurrently gives the correct result.
Thank you so much for the help (whoever you are), those are little subtle bugs that can drive you crazy. Now I know why our Prof. in Functional programming was making a big deal out of assignment statements! :) As for the domain between -pi/2 and pi/2 do I just keep adding/ subtracting until I get it within that domain? Thanks.
Chico wrote:
> On May 6, 7:41 pm, Nobody <nob...@nowhere.com> wrote: >> On Wed, 06 May 2009 14:33:14 -0700, Chico wrote: >>> Hello, I've been trying to figure out the way the CORDIC algorithm >>> works and I am not sure if this is the right place to post this, so >>> feel free to point me in the right direction to get help. >>> I usually code a quick prototype in python for verification and >>> understanding of an algorithm I need to implement. For those of you >>> who don't know Python but know C, you should be able to follow it >>> since python has a similar syntax to C. I've been struggling with this >>> CORDIC algorithm, seems straight forward from the wikipedia page, but >>> somehow I messing something up. >>> x = (x - (d*(2.0**(-i))*y)) >>> y = ((d*(2.0**(-i))*x) + y) >> You're calculating y[i+1] from x[i+1], when it should be x[i]. >> >> Rewriting the above as: >> >> (x,y) = (x - (d*(2.0**(-i))*y), (d*(2.0**(-i))*x) + y) >> >> to update x and y concurrently gives the correct result. > > Thank you so much for the help (whoever you are), those are little > subtle bugs that can drive you crazy. Now I know why our Prof. in > Functional programming was making a big deal out of assignment > statements! :) > > As for the domain between -pi/2 and pi/2 do I just keep adding/ > subtracting until I get it within that domain? Thanks.
Usually with embedded trig algorithms you just implement a single quadrant (0 to pi/2), but you can use a half-circle if you want. And yes, it's mostly a matter of adding/subtracting. But for some quadrants you might need to negate the angle, or the result - I'll let you figure out the details yourself. This is the first time I've heard someone describe Python syntax as being "similar to C" - it's about as different as you could get while still remaining an imperative programming language. But a well-written Python program should be fairly easily understood by C programmers. When posting sample code in a newsgroup (or other public forum), it's a good idea to cut out unnecessarily or old code - keep it to a minimum. Including the nested function K_vals(), for example, serves only to confuse.
On Wed, 06 May 2009 19:45:04 -0700, Chico wrote:

> As for the domain between -pi/2 and pi/2 do I just keep adding/ > subtracting until I get it within that domain? Thanks.
For any trigonometric function f, start by using the 2*pi periodicity, f(x+2*n*pi) = f(x) to coerce the value into the range 0..2*pi, e.g. x = x % 2*pi or x = x - 2*pi * floor(x / 2*pi) (but check how % handles negative values). Then, for values outside of the range 0..pi/2, use the identities: sin(x) = sin(pi - x) for pi/2 <= x <= pi sin(x) = -sin(x - pi) for pi <= x <= 3*pi/2 sin(x) = -sin(2*pi - x) for 3*pi/2 <= x <= 2*pi cos(x) = -cos(pi - x) for pi/2 <= x <= pi cos(x) = -cos(x - pi) for pi <= x <= 3*pi/2 cos(x) = cos(2*pi - x) for 3*pi/2 <= x <= 2*pi tan(x) = -tan(pi - x) for pi/2 <= x <= pi tan(x) = tan(x - pi) for pi <= x <= 3*pi/2 tan(x) = -tan(2*pi - x) for 3*pi/2 <= x <= 2*pi This ensures that the argument is always in the range 0..pi/2. IOW: Quadrant pre post 0..pi/2 n/a n/a pi/2..pi x = pi - x cos = -cos, tan = -tan pi..3*pi/2 x = x - pi sin = -sin, cos = -cos 3*pi/2..2*pi x = 2*pi - x sin = -sin, tan = -tan
In article <pan.2009.05.07.12.31.53.141000@nowhere.com>, 
nobody@nowhere.com says...
> On Wed, 06 May 2009 19:45:04 -0700, Chico wrote: > > > As for the domain between -pi/2 and pi/2 do I just keep adding/ > > subtracting until I get it within that domain? Thanks. > > For any trigonometric function f, start by using the 2*pi periodicity, > f(x+2*n*pi) = f(x) to coerce the value into the range 0..2*pi, e.g. > x = x % 2*pi or x = x - 2*pi * floor(x / 2*pi) (but check how % handles > negative values). > > Then, for values outside of the range 0..pi/2, use the identities: > > sin(x) = sin(pi - x) for pi/2 <= x <= pi > sin(x) = -sin(x - pi) for pi <= x <= 3*pi/2 > sin(x) = -sin(2*pi - x) for 3*pi/2 <= x <= 2*pi > > cos(x) = -cos(pi - x) for pi/2 <= x <= pi > cos(x) = -cos(x - pi) for pi <= x <= 3*pi/2 > cos(x) = cos(2*pi - x) for 3*pi/2 <= x <= 2*pi > > tan(x) = -tan(pi - x) for pi/2 <= x <= pi > tan(x) = tan(x - pi) for pi <= x <= 3*pi/2 > tan(x) = -tan(2*pi - x) for 3*pi/2 <= x <= 2*pi > > This ensures that the argument is always in the range 0..pi/2. > > IOW: > > Quadrant pre post > 0..pi/2 n/a n/a > pi/2..pi x = pi - x cos = -cos, tan = -tan > pi..3*pi/2 x = x - pi sin = -sin, cos = -cos > 3*pi/2..2*pi x = 2*pi - x sin = -sin, tan = -tan
Also known as the ASTC rule All, Sin, Tan, Cos For the quadrants from 0 (rad/degree) that results are positive for. -- Paul Carpenter | paul@pcserviceselectronics.co.uk <http://www.pcserviceselectronics.co.uk/> PC Services <http://www.pcserviceselectronics.co.uk/fonts/> Timing Diagram Font <http://www.gnuh8.org.uk/> GNU H8 - compiler & Renesas H8/H8S/H8 Tiny <http://www.badweb.org.uk/> For those web sites you hate
On May 7, 1:31=A0pm, Nobody <nob...@nowhere.com> wrote:
> On Wed, 06 May 2009 19:45:04 -0700, Chico wrote: > > As for the domain between -pi/2 and pi/2 do I just keep adding/ > > subtracting until I get it within that domain? Thanks. > > For any trigonometric function f, start by using the 2*pi periodicity, > f(x+2*n*pi) =3D f(x) to coerce the value into the range 0..2*pi, e.g. > x =3D x % 2*pi or x =3D x - 2*pi * floor(x / 2*pi) (but check how % handl=
es
> negative values). > > Then, for values outside of the range 0..pi/2, use the identities: > > =A0 =A0 =A0 =A0 sin(x) =3D sin(pi - x) =A0 =A0for pi/2 <=3D x <=3D pi > =A0 =A0 =A0 =A0 sin(x) =3D -sin(x - pi) =A0 for pi <=3D x <=3D 3*pi/2 > =A0 =A0 =A0 =A0 sin(x) =3D -sin(2*pi - x) for 3*pi/2 <=3D x <=3D 2*pi > > =A0 =A0 =A0 =A0 cos(x) =3D -cos(pi - x) =A0 for pi/2 <=3D x <=3D pi > =A0 =A0 =A0 =A0 cos(x) =3D -cos(x - pi) =A0 for pi <=3D x <=3D 3*pi/2 > =A0 =A0 =A0 =A0 cos(x) =3D cos(2*pi - x) =A0for 3*pi/2 <=3D x <=3D 2*pi > > =A0 =A0 =A0 =A0 tan(x) =3D -tan(pi - x) =A0 for pi/2 <=3D x <=3D pi > =A0 =A0 =A0 =A0 tan(x) =3D tan(x - pi) =A0 =A0for pi <=3D x <=3D 3*pi/2 > =A0 =A0 =A0 =A0 tan(x) =3D -tan(2*pi - x) for 3*pi/2 <=3D x <=3D 2*pi > > This ensures that the argument is always in the range 0..pi/2. > > IOW: > > =A0 =A0 =A0 =A0 Quadrant =A0 =A0 =A0 =A0pre =A0 =A0 =A0 =A0 =A0 =A0 post > =A0 =A0 =A0 =A0 0..pi/2 =A0 =A0 =A0 =A0 n/a =A0 =A0 =A0 =A0 =A0 =A0 n/a > =A0 =A0 =A0 =A0 pi/2..pi =A0 =A0 =A0 =A0x =3D pi - x =A0 =A0 =A0cos =3D -=
cos, tan =3D -tan
> =A0 =A0 =A0 =A0 pi..3*pi/2 =A0 =A0 =A0x =3D x - pi =A0 =A0 =A0sin =3D -si=
n, cos =3D -cos
> =A0 =A0 =A0 =A0 3*pi/2..2*pi =A0 =A0x =3D 2*pi - x =A0 =A0sin =3D -sin, t=
an =3D -tan At school we learned that the trig functions that are positive in each quadrant are: Quadrant +ve Mnemonic 0..pi/2 sin, cos, tan All pi/2..pi sin Stations pi..3*pi/2 tan To 3*pi/2..2*pi cos Crewe Might mean something if you are English ;->
Chico wrote:
> As for the domain between -pi/2 and pi/2 do I just keep adding/ > subtracting until I get it within that domain? Thanks.
If you represent angle such that 0 to pi maps to 0 to 2^n then getting within the domain only requires a bit-mask (no iteration). Bob
s0lstice wrote:
> Nobody <nob...@nowhere.com> wrote: >
... snip ...
> >> Quadrant pre post >> 0..pi/2 n/a n/a >> pi/2..pi x = pi - x cos = -cos, tan = -tan >> pi..3*pi/2 x = x - pi sin = -sin, cos = -cos >> 3*pi/2..2*pi x = 2*pi - x sin = -sin, tan = -tan > > At school we learned that the trig functions that are positive > in each quadrant are: > > Quadrant +ve Mnemonic > 0..pi/2 sin, cos, tan All > pi/2..pi sin Stations > pi..3*pi/2 tan To > 3*pi/2..2*pi cos Crewe
Much simpler, and not language dependent, is just think of the signs applied to the x and y axis when you draw a rt triangle in it. x is positive with displacement to right, y with displacement up. Now all you need to remember is what sin, cos, tan mean, i.e. the ratio of which sides. -- [mail]: Chuck F (cbfalconer at maineline dot net) [page]: <http://cbfalconer.home.att.net> Try the download section.
On Thu, 07 May 2009 16:51:10 -0400, CBFalconer <cbfalconer@yahoo.com>
wrote:

><snip> >Much simpler, and not language dependent, is just think of the >signs applied to the x and y axis when you draw a rt triangle in >it. x is positive with displacement to right, y with displacement >up. Now all you need to remember is what sin, cos, tan mean, i.e. >the ratio of which sides.
Your basic "unit circle" definitions. Sometimes, I wonder if it is even taught, anymore. Few younger folks in today's college calculus classes seem to even know about it -- and my latest experience with them was this last Winter term. Yet it's what I always go back to if I want to gain some insight on a new question (or one I have long since forgotten about.) You can't go wrong with that and Pythagorean. You also get all your half-angle and double-angle equivalences, etc. It's all there to be had with merely a gander. Also, there are symmetries to exploit beyond the 180 or 90 degree ones that I've seen mentioned in this thread. At a minimum, octants (or 45 degrees) are your basic brain-dead demarcation for calculation. More is still possible. For those whose approximation exposure has been sadly limited only to the very easily understood Taylor's theorem, then if for no other reason than the fact that staying closer in range to where the approximation was taken improves retained precision per cpu cycle consumed. Jon