Speeding up the CORDIC algorithm with a DSP

Digital signal processors (DSPs) crunch the numbers for applications that require fast analog-to-digital-to-analog conversion, such as software-defined radio and radar. If you weren't using a DSP, you might use the CORDIC algorithm to perform similar calculations on your IC. The strength of the CORDIC algorithm is its ability to solve vector rotation without using a multiplier.

Despite the CORDIC's well-documented properties, you won't often find it implemented on a DSP because the CORDIC was conceived 49 years ago when the cost of multiplier hardware was prohibitive. (DSPs are typically equipped with multipliers.) But what happens when you mix the two? When a CORDIC algorithm is implemented on a DSP processor, can the multipliers improve CORDIC's performance, or can they be left idle?

We're going to answer this question positively by proposing a novel way to map the CORDIC algorithm onto DSP hardware. This new method uses the DSP's MAC (multiply/accumulate) units and eliminates the conditionals, thus preserving the machine pipelining.

CORDIC, an acronym for COordinate Rotation DIgital Computer, is a class of shift-add algorithms that rotate a vector in a plane. CORDIC has become a commonly used method in memory- and CPU-constrained embedded systems because it's a simple and efficient way to calculate the hyperbolic and trigonometric functions found in every scientific calculator.

Most often the algorithm is used when no hardware multiplier is available, such as in microcontrollers and FPGAs, because the only operations it requires are addition, subtraction, bitshift, and table lookup.

CORDIC is widely used due to its simplicity and its property of relatively fast convergence. It has many applications, including computing trigonometric functions and converting Cartesian coordinates to polar coordinates (and vice versa).

Basically the CORDIC algorithm chooses special angles of rotation such that it can perform the rotation operations by simple shifts and additions, rather than the multiply functions that are required in the general case. Thus, design teams can use the CORDIC algorithm instead of hardware multipliers, which require a higher gate count and cost more to build.

Before we get into the details of how to map the CORDIC algorithm onto a DSP engine, we'll briefly review the CORDIC algorithm that calculates the magnitude and angle of a vector from its Cartesian coordinates. We'll then describe the implementation of the algorithm on a DSP processor, which uses addition and shift operations only.

First, however, let's review the traditional approach to implementing CORDIC.

Let v be a vector with Cartesian coordinates (x, y ). To simplify the description, we consider the right half plane of a unit circle only, in other words, we assume 1>x >0, 1>y .

The objective is to find the magnitude

and angle
.

If we can somehow rotate vector v = (x , y ) to v e = (x e , y e ) such that y e =0, the magnitude r will be the x -coordinate x e and the angle φ will be the rotated angle. This rotation is actually achieved by a number of successive rotations (called subrotations ). For each subrotation, the angle of rotation is properly chosen such that:

• The computation can be accomplished by addition and shift operations only (avoiding the use of multiplication).

• The set of subrotations will drive vector v to the x -axis, in other words, makingcoordinate y equal to 0; this can be guaranteed if the current rotation is half of the previous one.

Mathematically, the CORDIC algorithm can be described as follows:

Initially let:

x 0 =x

y 0 =y

φ0 =0

The first operation will rotate v0 = (x 0 , y 0 ) by α0 =45° to get v1 = (x 1 , y 1 ).

x 1 = x 0 cos(α0 ) – d 0 y 0 sin(α0 )        (1.3)

y 1 = y 0 cos(α0 ) + d 0 x 0 sin(α0 )        (1.4)

φ1 = φ0d 0 α0         (1.5)

where d 0 is the direction of rotation:

d 0 =1 if y 0 <0        (1.6)

d 0 =-1 if y 0 ≥0         (1.7)

The Equations 1.3 and 1.4 are the same as:

x 1 = cos(α0 )[x 0d 0 y 0 tan(α0 )]        (1.8)

y 1 = cos(α0 )[y 0 + d 0 x 0 tan(α0 )]        (1.9)

In general, the i -th rotation is:

x i +1 = cos(αi )[x i d i y i tan(αi )        (1.10)

y i +1 = cos(αi )[y i + d i x i tan(αi )]        (1.11)

φi +1 = φi d i αi                          (1.12)

where d i is the direction of rotation:

d i =1 if y i <0        (to rotate upwards)        (1.13)

d i =-1 if y i ≥0        (to rotate downwards)        (1.14)

for i =0,1,…

If chose ai such that tan(ai )= 2-i , the Equations 1.10 through 1.12 become:

x i +1 = K i [x i d i y i 2-i ]        (1.15)

y i +1 = K i [y i + d i x i 2-i ]        (1.16)

φi +1 = φi d i αi                           (1.17)

where:

K i = cos(arctan(2-i ))        (1.18)

αi = arctan(2-i )        (1.19)

Since the effect of constant K i can be removed later (see Equation 1.23), the CORDIC algorithm in Equations 1.15 through 1.17 can be implemented as:

x i +1 = x i d i y i 2-i         (1.20)

y i +1 = y i + d i x i 2-i         (1.21)

φi +1 = φi d i αi         (1.22)

where d i and αi are defined as in Equations 1.13, 1.14, and 1.19.

Suppose that, after N rotations, the desired precision is achieved, then the magnitude r and angle φ can then be found approximately by:

        (1.23)

        (1.24)

where K = K N -1 • …•K 0 = cos(arctan(2-( N -1) ))•… •cos(arctan(2-0 )), which can be precomputed precisely.

A typical pseudo code for CORDIC in Equations 1.20 through 1.24 is:

If (y i <0)        (1.25)

    x i +1 = x i y i '        (1.26)

    y i +1 = y i + x i '        (1.27)

    φi +1 = φi – αi         (1.28)

If (y i ≥ 0)        (1.29)

    x i +1 = x i + y i '        (1.30)

    y i +1 = y i x i '        (1.31)

    φi +1 = φi + αi         (1.32)

where:

x i '= x i 2-i = x i >>i         (1.33)

y i '= y i 2-i = y i >>i         (1.34)

i =0,1,…

the notation “z >>i ” means that the bits in z will be shifted down by i bits (the shift in this document is always an arithmetic shift, not a logical shift).

Figure 1 shows an illustration of the CORDIC algorithm. The pseudo code in Equations 1.25 through 1.32 is simple and requires only addition and shift operations. Since it doesn't require the multiply function, the hardware implementation is also simple, and for this reason the CORDIC algorithm is widely implemented in ASICs (application-specific integrated circuits) and FPGAs (field-programmable gate arrays).

View the full-size image

New CORDIC implementation for modern architectures
It's been a widely held belief that DSP processors can't improve the performance of the CORDIC algorithm. At first glance, this seems plausible, since the major CORDIC advantage in most implementations is its avoidance of multiply functions. However, if a multiplier is available, can it improve the performance of the CORDIC algorithm? We propose a reformulation of CORDIC that takes advantage of DSPs' hardware architecture .

There is another significant consideration when implementing CORDIC on a DSP architectures–that is the pipeline. To increase the performance, modern processors are usually pipelined often with as many as 10 stages. This pipelined design allows the processor to perform a task in a time-sliced fashion. A pipeline is an efficient approach as long as the code doesn't break the instruction flow with branches. The traditional CORDIC algorithm requires a conditional execution (see the pseudo code in Equations 1.25 through 1.32) that will typically break the pipeline and result in stalls. Hence, the performance of the CORDIC algorithm may be slow in a modern pipelined architecture.

Our new approach gets around this difficulty by mostly using multiplication units in a DSP processor; it also eliminates the need for conditional executions. We begin by reformulating the algorithm such that it's more suitable for using multipliers.

The iteration formula for x i +1 , y i +1 , and φi +1 (see Equations 1.25 to 1.32) may be written as:

x i +1 = φ i (x i y i ')+(1- f i ) (x i + y i ')        (2.1)

y i +1 = f i (y i + x i ')+(1- f i ) (y i x i ')        (2.2)

fi +1 = φ i i – ai )+(1- f i ) (φi + αi )        (2.3)

where:

f i =1 if y i <0        (2.4)

f i =0 if y i ≥0        (2.5)

x i '= x i 2-i =x i >>i        (2.6)

y i '=y i 2-i =y i >>i        (2.7)

The iteration in Equation 2.1 may be written as:

x i +1 = f i (– 2y i ')+ (x i + y i ')

  = 2[-f i y i '+ x i /2+ y i '/2]

  = 2[(-f i + 1/ 2)y i '+ x i /2]        (2.8)

  = 2[(-f i + 1/ 2)'y i + x i /2]

where:

(-f i + 1/ 2)' = (-f i + 1/ 2)2-i = (-f i + 1/ 2)>>i         (2.9)

Similarly, the iteration in Equations 2.2 and 2.3 may be written as:

y i +1 = 2[-(-f i + 1/ 2)'x i + y i /2]        (2.10)

φi +1 = 2i [(-f i + 1/ 2)αi i /2]        (2.11)

Collecting the above together, we have the algorithm to implement:

x i +1 = 2[(-f i + 1/ 2)'y i + x i /2]        (2.12)

y i +1 = 2[-(-f i + 1/ 2)'x i + y i /2]        (2.13)

φi +1 = 2i [(-f i + 1/ 2)α i i /2]        (2.14)

where:

f i =1 if y i <0        (2.15)

f i =0 if y i ≥ 0         (2.16)

(-f i + 1/ 2)' = (-f i + 1/ 2)2 i = (-f i + 1/ 2)>>i        (2.17)

The advantages of Equations 2.12 through 2.17 over Equations 1.25 through 1.32 are as follows:

1. The execution of Equations 2.12 through 2.17 is unconditional; hence, it gets around the difficulty of a broken pipeline, which was the case when implementing Equations 1.25 through 1.32.

2. The shift operation in Equations 2.12 through 2.17 will be done only once (to get the modified flag (-f i + 1/ 2)' ) while in Equations 1.25 through 1.32, two shift operations are needed (to get x i ' and y i '), it therefore saves one operation.

3. The flag has been chosen to have two possible values, 1/ 2 or –1/ 2. These values in Equation 1.15 fractional format are represented in hexadecimal notation as 0x4000 and 0xC000 respectively. Because of the particular values selection, the flags constants do not lose precision while shifting them down.

4. The new CORDIC formulation does not shift down the x i and y i coordinate values during the iteration process. Thus, no precision is lost on the original (x , y ) coordinate values.

5. The product of the flag, (-f i +1/ 2)', by x i or y i coordinate values is stored in a 40-bit accumulator. See Equations 2.8 and 2.10.

6. As a result of improvements 3, 4, and 5, the new CORDIC formulation achieves a higher precision of about 0.5 bit.

The formula in Equations 2.12 through 2.17 is particularly suitable to be implemented on a pipelined DSP architecture. When implemented on a Blackfin BF533 DSP processor from Analog Devices, each iteration took four cycles, while the traditional implementation in Equations 1.25 through 1.32 took seven cycles per iteration. For example, in a software radio application, we need to implement the CORDIC algorithm for two channels at 240 kHz. The number of iterations (subrotations) is 13. Our new method will consume 2×240,000x13x4 = 25 mips (million instructions per second) versus 2×240,000x13x7 = 44 mips.

This represents a saving of 43% in terms of mips. In addition, for the same number of iterations, the results of the new method in Equations 2.12 through 2.17 has, on average, a half bit more precision than the traditional method in Equations 1.25 through 1.32.

C & Assembly Code for CORDIC
In this section, we present three code samples:

Listing 1: C code implementing original CORDIC.

Listing 2: Blackfin assembly code implementing reformulated CORDIC.

Listing 3: Blackfin assembly code implementing original CORDIC.

The code in Listing 1 and the complete code in Listing 2 can be compiled and executed using ADI's VisualDSP++ tools. Note that the code in Listing 2 is an excerpt–the complete code is available at www.embedded.com/code/2008code.htm. The code assumes that |x |, |y |<1 in 1.15 format; in other words, x , y are in the range of [0x8000, 0x7fff] . The output phase has 32 bit in 1.31 format; in other words, -p is represented by 0x80000000 and p by 0x7fffffff . Assembly code may have some special requirement that may be seen in the comments.

View the full-size image

View the full-size image

For complete Listing 2 code, click here

View the full-size image

The code in Listing 3 is only a segment that shows how the original CORDIC may be implemented on ADI assembly. In Listing 3 , we show the code for one iteration only (total seven lines of code, costing seven cycles). ITER_NUM=i is a variable i =0,…,14 is the iteration number. It should be in a register, which may be loaded from memory. We leave it as it is for clarity. The register r2 = φi +1 , i =0,1,…,14 with r2=φ0 =0 initially (see 1.32). Besides the slow performance of the following implementation, it has less precision due to the line r1 = r0>>> ITER_NUM (v) , which shifts the x i and y i coordinates down, losing precision.

George Pan is a software engineer in the general purpose DSP Division of Analog Devices. He earned his degree in electrical engineering from the University of Connecticut. You may reach him at .

Fabian Lis is the Automotive Digital Radio engineering manager within DSPS automotive product line at ADI. Fabian holds a BSEE from Technion Israel Institute of Technology and an MSEE from Tel Aviv University. You may reach him at .

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.