Cordic, atan2, 32bit integerarithmetik, cortex M3

Gast #2114824 
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


1
//*****************************************************************************

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

3
//*****************************************************************************

4


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

6


7
  // Winkeltabelle: 90Grad,atan(1),atan(0.5),atan(0.25)...

8
  static const signed long i_atantab[] ={1073741824,536870912,316933405,

9
                167458907,85004756,42667331,21354465,10679838,5340245,

10
                2670163,1335087,667544,333772,166886,83443,41722,20861,

11
                10430,5215,2608,1304,652,326,163,81,41,20,10,5,3,1};//31

12


13
  // Eingangswerte begrenzen

14
  #define lmax (0x3FFFFFFF) 

15
  while ( x>lmax || (-x)>lmax || y>lmax || (-y)>lmax ) {

16
  x = x>>1 ; y = y>>1;

17
  }

18


19
  long sign;

20
  long k=0;

21
  long xh,yh;

22
  

23
  // initalisierung

24
  sign = (y>=0) ? -1 : 1;

25
  xh=-sign*y;  yh=+sign*x;

26
  x=xh;y=yh;

27
  *winkel=sign*i_atantab[0];

28


29
  // iteration maximum: i<31 

30
  for (int i=1;i<15;i++) {

31
    sign = (y>=0) ? -1 : 1;

32
    *winkel += sign*i_atantab[i];

33
    xh=x-sign*(y>>k);  

34
    yh=y+sign*(x>>k);

35
    x=xh;y=yh;

36
    k=k+1;

37
  }

38


39
}
Antwort Zitieren Bearbeiten Thread verschieben Thread sperren Langsamer Modus Anmeldepflicht aktivieren Anpinnen Thread löschen Thread mit anderem zusammenführen Melden Danke
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
Antwort Zitieren Bearbeiten Löschen Melden Danke
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
Antwort Zitieren Bearbeiten Löschen Melden Danke 
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:

1
//*****************************************************************************

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

3
//*****************************************************************************

4


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

6


7
  // Winkeltabelle: 90Grad,atan(1),atan(0.5),atan(0.25)...

8
  static const signed long i_atantab[] ={1073741824,536870912,316933405,

9
                167458907,85004756,42667331,21354465,10679838,5340245,

10
                2670163,1335087,667544,333772,166886,83443,41722,20861,

11
                10430,5215,2608,1304,652,326,163,81,41,20,10,5,3,1};//31

12


13
  // Eingangswerte begrenzen

14
  #define lmax (0x1FFFFFFF) 

15
  long sign;

16
  long k=0;

17
  long xh,yh,hy,winkel;

18
  int  nshift=0;  

19
  while ( x>lmax || (-x)>lmax || y>lmax || (-y)>lmax ) {

20
  x = x>>1 ; y = y>>1; nshift++;

21
  }

22
  

23
  // initalisierung

24
  sign = (y>=0) ? -1 : 1;

25
  xh=-sign*y;  yh=+sign*x;

26
  x=xh;y=yh;

27
  winkel=sign*i_atantab[0];

28


29
  // iteration maximum: i<31 

30
  for (int i=1;i<=15;i++) {

31
    sign = (y>=0) ? -1 : 1;

32
    winkel += sign*i_atantab[i];

33
    xh=x-sign*(y>>k);  

34
    yh=y+sign*(x>>k);

35
    x=xh;y=yh;

36
    k=k+1;

37
  }

38
  // auf x bleibt die hypotenuse stehen, muss aber noch mit

39
  // einem zur Iterationstiefe passenden Faktor scaliert werden.

40
  // siehe http://en.wikipedia.org/wiki/CORDIC

41
  // für 15 ist da 0,60725293538591 angegeben

42
  if (hypotenuse)

43
    { 

44
    #define SCALE_FAKTOR ((long long int)(0.60725293538591 * (1<<30)))

45
    hy=(long)(((long long int)x * SCALE_FAKTOR)>>30);

46
    hy<<=nshift;  

47
    *hypotenuse=hy;

48
    }

49
  return -winkel;

50
}


Habe ich das mit der Scalierungskonstanten richtig verstanden?
Antwort Zitieren Bearbeiten Löschen Melden Danke 
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

1
import java.util.Locale;

2


3
public class atan2_v2 {

4


5
  static int i_atantab[] = {1073741824,536870912,316933405,

6
    167458907,85004756,42667331,21354465,10679838,5340245,

7
    2670163,1335087,667544,333772,166886,83443,41722};

8


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

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

11


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

13


14
    // Generierung der Tabelle: TODO: aus dem c-code loeschen

15
    //    System.out.printf("static const signed long i_atantab[] ={%d",w90);

16
    //    double arg=1;

17
    //    for (int i=1;i<31;i++) {

18
    //      int atan = (int) Math.round( Math.atan(arg) * w180 / Math.PI );

19
    //      System.out.printf(",%d",atan);

20
    //      arg=arg/2;

21
    //    }

22
    //    System.out.printf("};\n");

23


24


25
    int sign;

26
    int k=0;

27
    int xh,yh;

28
    int winkel;

29
    int lmax = (0x1FFFFFFF); 

30


31
    // die Eingangswerte begrenzen

32
    int  nshift=0;  

33
    while ( x>lmax || (-x)>lmax || y>lmax || (-y)>lmax ) {

34
      x = x>>1 ; y = y>>1; nshift++;

35
    }

36


37
    // initalisierung

38
    sign = (y>=0) ? -1 : 1;

39
    xh =- sign*y;  yh =+ sign*x;

40
    x=xh; y=yh;

41
    winkel = sign*i_atantab[0];

42


43
    // iteration maximum: i<31 

44
    for (int i=1;i<14;i++) {

45
      sign = (y>=0) ? -1 : 1;

46
      winkel += sign*i_atantab[i];

47
      //  System.out.printf(Locale.ENGLISH,

48
      //      "i_atan2: %d: x=%d ; y=%d ; sig=%d ; k=%d ; res=%f\n",

49
      //      i-1, x,y,sign,k,(double)(res*180.0/w180));

50
      xh=x-sign*(y>>k);  

51
      yh=y+sign*(x>>k);

52
      x=xh;y=yh;

53
      k=k+1;

54
    }

55


56


57
    // auf x bleibt die hypotenuse stehen, muss aber noch mit

58
    // einem zur Iterationstiefe passenden Faktor scaliert werden.

59
    // siehe http://en.wikipedia.org/wiki/CORDIC

60
    // für 15 ist da 0,60725293538591 angegeben

61
    if (hypotenuse>0)  { 

62
      int SCALE_FAKTOR =  (int)  ( 0.60725293538591 * (1<<30));

63
      hypotenuse = (int) ( ((long) x * SCALE_FAKTOR)>>30 );

64
      hypotenuse <<= nshift;  

65
      System.out.println("hypotenuse=" + hypotenuse);

66
    }

67


68
    return -winkel;

69
  }

70


71


72
  /** ***************************************************

73
   * @param args

74
   */

75
  public static void main(String[] args) {

76
    int x=2*10000;

77
    int y=7*10000;

78
    double res1 = Math.atan2(y, x);

79


80
    int res2 = i_atan2(1,x,y);

81
    System.out.printf(Locale.ENGLISH,

82
        "x=%d ; y=%d : w1=%3.5f : w2=%3.5f\n", 

83
        x,y,res1*180.0/Math.PI,(double)res2*180.0/w180);

84


85
  }

86


87
}
Antwort Zitieren Bearbeiten Löschen Melden Danke
Thread beobachten

Antwort schreiben

Eine eigene Frage? Dieser Thread ist über 15 Jahre alt. Stell sie in einem eigenen Thread – dort lesen mehr Leute mit, und sie wird öfter und schneller beantwortet.

Neuen Thread starten

Eine Lösung, Ergänzung oder Rückmeldung zu diesem Thread? Dann antworte hier.

Formatierung
  • [c]C-Code[/c]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] – Link zu einem Artikel
  • Verweis auf einen anderen Beitrag: Rechtsklick auf dessen Beitrags-ID (#8078628), „Adresse kopieren“, in den Text einfügen

Mehr zur Formatierung

Längeren Sourcecode als Anhang, nicht in den Text. JPG nur für Fotos und Scans; Zeichnungen und Screenshots als PNG oder GIF. Siehe Bildformate.

Dein Name und deine E-Mail-Adresse

Schon registriert? Anmelden

Wird nicht angezeigt. Wir schicken dir einen Link; sobald du damit bestätigst, erscheint dein Beitrag.

Mit dem Abschicken legst du einen kostenlosen Account an und erkennst die Nutzungsbedingungen und die Datenschutzerklärung an.