DSP Main | Archives | Feedback

Take a look at any DSP algorithm, and you're likely to see sine and cosine functions embedded somewhere in the code. Some obvious examples include FFTs and numerically controlled oscillators. It's easy to implement the constituent trig functions in software using lookup tables or polynomial approximations, but what alternatives are available if a realtime problem demands a high-speed solution in hardware?

In such cases, readers might be interested to know about a class of iterative algorithms (Ref 1) that implement many useful functions -- among them the sine and cosine -- but involve nothing more complicated than bit shifts and additions. This technique, called the CORDIC (COordinate Rotation DIgital Computer) algorithm, yields simple hardware implementations that fit nicely into an FPGA or ASIC.

The basis of this technique is the rotation matrix

Given a 2D vector x, A(θ)x rotates x by a positive angle θ. If the two components of x correspond to the real and imaginary parts of a complex number, then A(θ)x corresponds to a multiplication of the complex number by a factor ejθ. This operation is equivalent to a phase shift of a complex waveform by an angle θ, the basis of many DSP algorithms including FFTs, modulators and demodulators.

Now consider where the CORDIC technique comes into play. It evaluates A(θ) in terms of tan θ and restates the above equation as

This matrix multiplication is particularly simple to implement if you choose θ = tan^-12^-k :

Ignoring for the moment the gain factor , it's apparent that because the matrix consists only of powers of two that computing A(θ)x requires nothing more than bit shifts and adds. The CORDIC technique cascades these simple matrix multiplications as shown in the following equation to approximate a rotation by an arbitrary angle as

A(θk) = A(±θ1)A(±θ2)...A(±θn )x,

where the algorithm chooses the signs of the angles to render increasingly good approximations to A(θ)x. Note that the normalization factors

don't depend on the value of θ, only on the number of stages n. Thus you can precompute the product of them once or even ignore them if a fixed scale factor is acceptable.

Fig 1 --Twelve iterations of the CORDIC algorithm yield a worst-case error of approximately 0.05 deg.

A natural question to ask concerns the accuracy of this approximation. To investigate this aspect, the Matlab code in Listing 1 implements the CORDIC rotation algorithm. The code starts from an initial test vector x = [1,0]T, performs a CORDIC rotation by -180º, -179º, ... , +180º and then computes the residual errors, which appear in Fig 1. The gain factor and atan() computations don't depend on the desired angle of rotation, so a hardware implementation could use precomputed values. Note that the variable delta is nothing more than 2-k, and thus a hardware implementation could easily achieve delta*x or delta*y as a bit shift.

It now should be apparent why CORDIC is a handy method to have when performing phase shifts in simple digital hardware for high-speed DSP applications. A simple modification of the technique (left to the reader as an exercise) inverts the function so that when given a vector [u,v]T, the algorithm returns the magnitude and phase tan-1 v/u.

Reference

Volder, JE, "The CORDIC trigonometric computing technique," IRE Transactions on Electronic Computers, 1959, pgs 330-334.

Mark Sullivan (dalek@radix.net) is Chief Scientist at SkyBitz Inc (Herndon, VA), a developer of tracking and communications services based on the GLS (Global Locating Systems) technology it invented. Mark received a PhD in Information Technology from George Mason Univ.

This article originally appeared in Personal Engineering & Instrumentation News, November 1998, pgs 56-58. Reprinted with permission of PEC Inc; all rights reserved.

%  Listing 1 -- test of the CORDIC algorithm for phase rotation

clear; 

n=12;                % number of iterations

for i=-180:180

   phi(i+181)=i;

   theta=phi(i+181)*pi/180;

   x=1;

   y=0;      %  start with a vector at zero

   % degrees the first iteration performs a special rotation by +/- 90 deg

   if (theta>pi/2)

      theta=theta-pi/2;

      newx=-y;

      newy=x;

   elseif (theta<-pi/2)

      theta=theta+pi/2;

      newx=y;

      newy=-x;

   else

      newx=x;

      newy=y;

   end

   x=newx;

   y=newy;

   delta=1;

   gain=1;

   % now perform the remaining iterations of CORDIC

   for k=2:n

      if(theta>0)

         theta=theta-atan(delta);

         newx=x-y*delta;

         newy=y+x*delta;

      else

         theta=theta+atan(delta);

         newx=x+y*delta;

         newy=y-x*delta;

      end

      x=newx;

      y=newy;

      gain=gain*sqrt(1+delta*delta);

      % account for growth of vector

      delta=delta/2;

      end

   x=x/gain;

   y=y/gain;

   error(i+181)=atan2(y,x)*180/pi-phi(i+181);

   if(error(i+181)<-180)

      error(i+181)=error(i+181)+360;

   elseif

      (error(i+181)>180)

      error(i+181)=error(i+181)-360;

   end

end

plot(phi,error);

		Click here to get your listing up.