Cordic, atan2, 32bit integerarithmetik, cortex M3

[OR (Gast)]

Hallo,

ich habe für 32bit den Cordic-algorithmus implementiert. Wer es mag, 
kann bis zu letzten Stelle rechnen, ich habe nach den ersten 15bit schon 
abgebrochen. Das Ergebnis ist ein 32bit Winkel (-180Grad=-2^31, 
180Grad=+2^31).

Falls jemand Bedarf dafür hat, einfach kopieren und einbauen.

Ich denke, der Code sollte schon recht rechenzeitoptimal sein. 
Vielleicht fallen Euch ja noch ein paar Optimierungen ein.

Viel Spaß damit
OR

//*****************************************************************************
// Funktion zum Errechnen des atan2. Cordic Algorithmus (OR2011)
//*****************************************************************************

void atan2_c(long *winkel , long x , long y) {

  // Winkeltabelle: 90Grad,atan(1),atan(0.5),atan(0.25)...
  static const signed long i_atantab[] ={1073741824,536870912,316933405,
                167458907,85004756,42667331,21354465,10679838,5340245,
                2670163,1335087,667544,333772,166886,83443,41722,20861,
                10430,5215,2608,1304,652,326,163,81,41,20,10,5,3,1};//31

  // Eingangswerte begrenzen
  #define lmax (0x3FFFFFFF) 
  while ( x>lmax || (-x)>lmax || y>lmax || (-y)>lmax ) {
  x = x>>1 ; y = y>>1;
  }

  long sign;
  long k=0;
  long xh,yh;
  
  // initalisierung
  sign = (y>=0) ? -1 : 1;
  xh=-sign*y;  yh=+sign*x;
  x=xh;y=yh;
  *winkel=sign*i_atantab[0];

  // iteration maximum: i<31 
  for (int i=1;i<15;i++) {
    sign = (y>=0) ? -1 : 1;
    *winkel += sign*i_atantab[i];
    xh=x-sign*(y>>k);  
    yh=y+sign*(x>>k);
    x=xh;y=yh;
    k=k+1;
  }

}

[Serdar U. (serdar)]

Hallo,

habe versucht das obige zu verwenden ...
Die arctan2(1;1) sollte  in Grad pi/4 - Radiant (0,785398) liefern.

Jedoch erhält man mit der obigen Funktionen einen ganz anderen Wert...

Ist der Code wirklich korrekt ??

[OR (Gast)]

Hallo,

die Funktion ist auf Festkommaarithmetik ausgelegt. Du solltest eher 
etwas wie artan2(&winkel,10000,10000) versuchen.

Grüße
OR

[Erich (Gast)]

Hallo "OR",
die Funktion realisiert den CORDIC Algorithmus.
Prinzipiell erscheint sie mir in Ordnung zu sein,
jedoch sind die Werte in der Tabelle nicht alle richtig.
Ich habe es mit Excel nachgerechnet.
Demnach ist der erste Wert ok, die folgenden 10 jedoch nicht ganz.
Ab dem 12. Wert passt es wieder.

Meine berechneten Werte sind

1.073.741.824
  633.866.811
  334.917.815
  170.009.512
   85.334.662
   42.708.931
   21.359.677
   10.680.490
    5.340.327
    2.670.173
    1.335.088
      667.544
      333.772
      166.886
      83.443
      41.722
      20.861
      10.430
       5.215
       2.608
       1.304
         652
         326
         163
          81
          41
          20
          10
           5
           3
           1

[Vladimir Baykov (Gast)]

You can improve the accuracy of arctan  if use Meggitt algorithm


Dr. Meggitt was the first who suggested the algorithms similar to 
CORDIC,
but using normalised constants and normalised variables reducible to 
zero.
It allowed him to improve significantly the accuracy of the iterative 
al-
gorithms.
In addition relative errors of the fixed-point representation numbers in
Meggitt's algorithms in contrast with Volder's algorithms practicaly
don't depend on the word size, like in the floating-point number
representation.
In other words, Dr.Meggitt suggested not only "pseudodivision and 
pseudo-
multiplication processes" (as indicated in the title of his paper), but 
also
"pseudofloating- point number representation" for CORDIC method.

You can find the Meggitt*s algorithms and results of its errors in
http://baykov.de/CORDIC1972.htm

Vladimir Baykov

[Jürgen L. (jliegner)]

Hallo und vielen Dank, ich habe das mal probiert. Der Winkel ist vom 
Vorzeichen falsch. Ausserdem werden falsche Ergebnisse berechnet wenn x 
und y in die Nähe von 0x3fffffff kommen. Ausserdem wäre es noch 
wünschenswert gleich die Hypotenuse mit zu bekommen. Folgende Vorschläge 
zum Code deshalb von mir:

//*****************************************************************************
// Funktion zum Errechnen des atan2. Cordic Algorithmus (OR2011)
//*****************************************************************************

long atan2_c(long *hypotenuse, long x , long y) {

  // Winkeltabelle: 90Grad,atan(1),atan(0.5),atan(0.25)...
  static const signed long i_atantab[] ={1073741824,536870912,316933405,
                167458907,85004756,42667331,21354465,10679838,5340245,
                2670163,1335087,667544,333772,166886,83443,41722,20861,
                10430,5215,2608,1304,652,326,163,81,41,20,10,5,3,1};//31

  // Eingangswerte begrenzen
  #define lmax (0x1FFFFFFF) 
  long sign;
  long k=0;
  long xh,yh,hy,winkel;
  int  nshift=0;  
  while ( x>lmax || (-x)>lmax || y>lmax || (-y)>lmax ) {
  x = x>>1 ; y = y>>1; nshift++;
  }
  
  // initalisierung
  sign = (y>=0) ? -1 : 1;
  xh=-sign*y;  yh=+sign*x;
  x=xh;y=yh;
  winkel=sign*i_atantab[0];

  // iteration maximum: i<31 
  for (int i=1;i<=15;i++) {
    sign = (y>=0) ? -1 : 1;
    winkel += sign*i_atantab[i];
    xh=x-sign*(y>>k);  
    yh=y+sign*(x>>k);
    x=xh;y=yh;
    k=k+1;
  }
  // auf x bleibt die hypotenuse stehen, muss aber noch mit
  // einem zur Iterationstiefe passenden Faktor scaliert werden.
  // siehe http://en.wikipedia.org/wiki/CORDIC
  // für 15 ist da 0,60725293538591 angegeben
  if (hypotenuse)
    { 
    #define SCALE_FAKTOR ((long long int)(0.60725293538591 * (1<<30)))
    hy=(long)(((long long int)x * SCALE_FAKTOR)>>30);
    hy<<=nshift;  
    *hypotenuse=hy;
    }
  return -winkel;
}

Habe ich das mit der Scalierungskonstanten richtig verstanden?

[Oliver R. (oliverr)]

Hallo,

das Vorzeichen habe ich nie ausprobiert. Danke.
Die Idee mit der Hypothenuse ist gut. Zum Testen des codes habe ich das 
ganze eigentlich in Java geschrieben. Ich habe Deinen Code mit 
reingebaut.

Vielen Dank für Deine Mühe
Oliver

import java.util.Locale;

public class atan2_v2 {

  static int i_atantab[] = {1073741824,536870912,316933405,
    167458907,85004756,42667331,21354465,10679838,5340245,
    2670163,1335087,667544,333772,166886,83443,41722};

  static int w180 = (int) Math.pow(2,31);
  static int w90 = (int) Math.pow(2,30);

  static int i_atan2(int hypotenuse, int x,int y) {

    // Generierung der Tabelle: TODO: aus dem c-code loeschen
    //    System.out.printf("static const signed long i_atantab[] ={%d",w90);
    //    double arg=1;
    //    for (int i=1;i<31;i++) {
    //      int atan = (int) Math.round( Math.atan(arg) * w180 / Math.PI );
    //      System.out.printf(",%d",atan);
    //      arg=arg/2;
    //    }
    //    System.out.printf("};\n");


    int sign;
    int k=0;
    int xh,yh;
    int winkel;
    int lmax = (0x1FFFFFFF); 

    // die Eingangswerte begrenzen
    int  nshift=0;  
    while ( x>lmax || (-x)>lmax || y>lmax || (-y)>lmax ) {
      x = x>>1 ; y = y>>1; nshift++;
    }

    // initalisierung
    sign = (y>=0) ? -1 : 1;
    xh =- sign*y;  yh =+ sign*x;
    x=xh; y=yh;
    winkel = sign*i_atantab[0];

    // iteration maximum: i<31 
    for (int i=1;i<14;i++) {
      sign = (y>=0) ? -1 : 1;
      winkel += sign*i_atantab[i];
      //  System.out.printf(Locale.ENGLISH,
      //      "i_atan2: %d: x=%d ; y=%d ; sig=%d ; k=%d ; res=%f\n",
      //      i-1, x,y,sign,k,(double)(res*180.0/w180));
      xh=x-sign*(y>>k);  
      yh=y+sign*(x>>k);
      x=xh;y=yh;
      k=k+1;
    }


    // auf x bleibt die hypotenuse stehen, muss aber noch mit
    // einem zur Iterationstiefe passenden Faktor scaliert werden.
    // siehe http://en.wikipedia.org/wiki/CORDIC
    // für 15 ist da 0,60725293538591 angegeben
    if (hypotenuse>0)  { 
      int SCALE_FAKTOR =  (int)  ( 0.60725293538591 * (1<<30));
      hypotenuse = (int) ( ((long) x * SCALE_FAKTOR)>>30 );
      hypotenuse <<= nshift;  
      System.out.println("hypotenuse=" + hypotenuse);
    }

    return -winkel;
  }


  /** ***************************************************
   * @param args
   */
  public static void main(String[] args) {
    int x=2*10000;
    int y=7*10000;
    double res1 = Math.atan2(y, x);

    int res2 = i_atan2(1,x,y);
    System.out.printf(Locale.ENGLISH,
        "x=%d ; y=%d : w1=%3.5f : w2=%3.5f\n", 
        x,y,res1*180.0/Math.PI,(double)res2*180.0/w180);

  }

}