The C Users' Group Library 1994 August/wc-cdrom-cusersgrouplibrary-1994-08.iso/listings/v_11_09/cordic.cpp

C/C++ Source or Header | 1993-02-28 | 5KB | 191 lines

//CORDIC.CPP  the implementation of methods prototyped in CORDIC.H

#include "CORDIC.HPP"

#define TEST


                //define static members declared in the class declaration
int cordic::ArcTan[NBits];  //array of integer arcTan's
int cordic::xInit = 0;          //length of initial vector
int cordic::instance = 0;             //previous instance of CORDIC?
long cordic::CordicBase = 0;
long cordic::HalfBase = 0;
long cordic::Quad2Boundary = 0;
long cordic::Quad3Boundary = 0;

                //define a handy constant
const long sixtyfourK = 0x10000;


void cordic::cordic(void) //constructor:  initializes statics for use later
{
    if(!instance)          //only one instance
    {
        int i;                //index arcTan[]
        double f;             //calc initial x projection
        long powr;        //calc powers of 2 up to 2^(2*(NBits-1))

        CordicBase = 1 << NBits;      //(* 2^NBits)
        HalfBase = CordicBase >> 1;   //(/2)
        Quad2Boundary = CordicBase << 1;
        Quad3Boundary = CordicBase + Quad2Boundary;

        //a diminishing series of ArcTan's
        powr = 1;
        for(i = 0; i < NBits; i ++)
        {
            ArcTan[i] = int(atan(1.0/powr)/(M_PI/2)*CordicBase + 0.5);
            powr <<= 1;
        }

        //figure the starting x by compensating for elongation during the
        //expansion of the series.
        f = 1.0;
        powr = 1;
        for (i = 0;i < NBits;i++)
        {
            f = (f * (powr + 1)) / powr;
            powr <<= 2;                     //even powers of 2
        }
        f = 1.0/sqrt(f);
        xInit = int(CordicBase * f + 0.5);
    }// if(!instance)
    instance ++;
}//cordic::cordic()

cordic::~cordic()  //destructor--does nothing in this implementation
{
    if(!--instance){}
}

long cordic::cordicbase(void) { return CordicBase;}

//the following routine takes input only in positive cordic angle units
//within the first 64K(once around the circle)
void cordic::sinCos(unsigned long theta, int &sin, int &cos)    // input in cordic Units
{
    int quadrant;     //quadrant of input theta
    int z;            //theta in the first quadrant
    int i;            //index vector rotations
    int x, y, x1, y1; //projections before and after rotation

    //determine quadrant of input angle theta, translate to 1st quadrant,
    //recording original quadrant for adjusting sign at end of computation
    if(theta < CordicBase)
    {
        quadrant = Quad1;
        z = int(theta);
    }
    else if(theta < Quad2Boundary)
    {
        quadrant = Quad2;
        z = int(Quad2Boundary - theta);
    }
    else if(theta < Quad3Boundary)
    {
        quadrant = Quad3;
        z = int(theta - Quad2Boundary);
    }
    else
    {
        quadrant = Quad4;
        z = - int(theta);
    }

    //initialize the projections on the x and y axes
    x = xInit;
    y = 0;

    //negate z so target angle will be 0 while vector rotates in 1st quad.
    z = -z;

    //rotate nBits times.
    for(i = 0;i < NBits; i++)
    {
        if(z < 0)
        {
            //Counterclockwise rotation
            z += ArcTan[i];
            y1 = y + (x >> i);
            x1 = x - (y >> i);
        }
        else
        {
            //Clockwise rotation
            z -= ArcTan[i];
            y1 = y - (x >> i);
            x1 = x + (y >> i);
        }

        //put new projections into old for next iter.
        x = x1;
        y = y1;
    }// for i = 0 .. nBits

    //determine sign based on quadrant
    cos = (quadrant == Quad1 || quadrant == Quad4) ? x : -x;
    sin = (quadrant == Quad1 || quadrant == Quad2) ? y : -y;
}//cordic::sinCos(long, int, int)


//this function takes inputs in cordic angle units greater than 64K, and
//less than 0, up to the size of long
void cordic::sinCos(long theta, int &sin, int &cos) //input of any long
{
    theta %= sixtyfourK;
    if(theta < 0)
        theta += sixtyfourK;
    sinCos((unsigned long)theta, sin, cos); //1st ver. of sinCos does the work
}
//takes input in positive or negative radians
void cordic::sinCos(float theta, int &sin, int &cos)   // input in Radians
{
    long CordicTheta;
    CordicTheta = long((theta / (2 * M_PI)) * sixtyfourK);
    sinCos(CordicTheta, sin, cos);  //second version of sinCos does the work
}

//define predefined instance
cordic cord;


#ifdef TEST
#include <iostream.h>
int main()
{
    unsigned long thetau;
    long theta;
    float thetaf;
    int sine, cosine;
    cout << endl << "The unsigned long version: " << endl;
    for(thetau = 0;thetau <= sixtyfourK; thetau += 300)
    {
        cord.sinCos(thetau, sine, cosine);
        cout << thetau << "  "
                << sine << "  "
                << cosine << endl;
    }//for(thetau...)
    cout << endl << "The long version: " << endl;
    for(theta = - sixtyfourK / 2;theta <= sixtyfourK * 2 ; theta += 800)
    {
        cord.sinCos(theta, sine, cosine);
        cout << theta << "  "
                << sine << "  "
                << cosine << endl;
    }//for(theta...)
    cout << endl << "The float version: " << endl;
    thetaf = 0;
    for(theta = 0;theta <= 360 ; theta++)
    {
        cord.sinCos(thetaf, sine, cosine);
        cout << (thetaf * 180 / M_PI) << "  "
                << sine << "  "
                << cosine << endl;
        thetaf += (2 * M_PI)/360;
    }//for(theta...)  thetaf
    return 0;
}


#endif //TEST