Cordic, atan2, 32bit integerarithmetik, cortex M3

OR Gast #2114824
Lesenswert?
• ▲
▼ 
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
}
Beitrag melden Bearbeiten Thread verschieben Thread sperren Anmeldepflicht aktivieren Anpinnen Thread löschen Thread mit anderem zusammenführen Markierten Text zitieren Antwort Antwort mit Zitat
Erich Gast #2140187
Lesenswert?
• ▲
▼ 
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
Beitrag melden Bearbeiten Löschen Markierten Text zitieren Antwort Antwort mit Zitat
Vladimir Baykov Gast #2141305
Lesenswert?
• ▲
▼ 
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
Beitrag melden Bearbeiten Löschen Markierten Text zitieren Antwort Antwort mit Zitat
Jürgen #2577279
Lesenswert?
• ▲
▼ 
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?
Beitrag melden Bearbeiten Löschen Markierten Text zitieren Antwort Antwort mit Zitat
Oliver R. #2580405
Lesenswert?
• ▲
▼ 
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
}
Beitrag melden Bearbeiten Löschen Markierten Text zitieren Antwort Antwort mit Zitat
Thread beobachten

Antwort schreiben

Bitte melde dich an, um einen Beitrag zu schreiben.

Bestehender Account

oder

Mit Google-Account einloggen

Noch kein Account?

Die Registrierung ist kostenlos und dauert nur eine Minute.

Jetzt registrieren