CORDIC.

En septembre 1959, Jack E. Volder publia "The CORDIC Trigonometric Computing Technique" (IRE Transactions on Electronic Computers). Il y décrivait une méthode itérative permettant de calculer les fonctions trigonométriques sinus, cosinus et les 'amplitude et phase' d'un vecteur. La particularité de cette méthode est de ne pas nécessiter de multiplication autre que par 2 et est par conséquent appropriée pour certains composants électroniques (par exemple, un microcontôleur ne possédant de multiplication hardware).

Note: "CORDIC" est un acronyme de COordinate Rotation DIgital Computer.

La méthode est basée sur des 'rotations' successives dont la somme des angles de rotation approche de plus en plus l'angle considéré. Cette approche peut aussi être vue comme des produits de nombres complexes puisqu'on sait que: 
si x=r (cos t + i sin t) et x'=r' (cos t' + i sin t')
alors x*x'=r*r' (cos (t+t') + i sin (t+t'))

ou exprimé autrement: 

après la multiplication de x par x', 
la partie réelle de x:      r cos t devient rr' cos (t+t')
la partie imaginaire de x:  r sin t devient rr' sin (t+t')

Si r'=1, seule la phase/argument change et le module/amplitude reste inchangé.
x*x'=r (cos (t+t') + i sin (t+t'))
Pour rappel, la transformation de x=a+bi en x=r(cos t + i sin t) est décrite ici.

L'idée de J.E.Volder a été de multiplier le nombre complexe initial par des nombres 1+Ki avec K=1, 2-1, 2-2, ..., 2-n, ...
où K=tg t et donc l'angle de rotation dont on parlait plus haut t=arctg K.
Malheureusement, le nombre 1+Ki n'a pas un module unitaire et par conséquent, il apparaît un gain à chaque multiplication.

Voici un tableau avec les différentes valeurs:
n K = 2-n x = 1 + Ki Argument de x
en degrés
= arctg(K)		 Module de x
=SQRT(1+K²) 		Gain
=module de x
*gain précédent
0 1.0 1 + 1.0 i 45.00000 1.41421356 1.414213562
1 0.5 1 + 0.5 i 26.56505 1.11803399 1.581138830
2 0.25 1 + 0.25 i 14.03624 1.03077641 1.629800601
3 0.125 1 + 0.125 i  7.12502 1.00778222 1.642484066
4 0.0625 1 + 0.0625 i  3.57633 1.00195122 1.645688916
5 0.03125 1 + 0.031250 i  1.78991 1.00048816 1.646492279
6 0.015625 1 + 0.015625 i  0.89517 1.00012206 1.646693254
7 0.007813 1 + 0.007813 i  0.44761 1.00003052 1.646743507
... ... ... ... ... ... 

Si x'=1+Ki, après multiplication de x=r (cos t + i sin t) par x',
la partie réelle de x:      r cos t devient r G cos (t+arctg K)
la partie imaginaire de x:  r sin t devient r G sin (t+arctg K)
où le gain G = module de (1+Ki) = SQRT(1+K²)
mais nous avons aussi (cf (a+bi)*(a'+b'i)=(aa'-bb')+i(ab'+a'b)):
la partie réelle de x:      r cos t devient r (cos t - K sin t)
la partie imaginaire de x:  r sin t devient r (K cos t + sin t)
Supposons devoir calculer sin 53° et cos 53°. Commençons avec un angle de 0° que nous allons progressivement faire approcher de 53°.
x=cos 0° + i sin 0° = 1+i 0
Multiplions par:
1+Ki avec K=1 (correspondant à un angle de 45°). Pourquoi 1+Ki et pas 1-Ki? Il faut augmenter l'angle actuel (0°) jusqu'à arriver à 53°.
Nous obtenons:
un angle de 0°+45°=45° (manquent 8°)
une partie réelle=1.414213562 cos 45° mais aussi = 1-K*0 = 1
une partie imaginaire=1.414213562 sin 45° mais aussi = K*1+0 = 1
         cf x=(a+bi)*(a'+b'i)=(aa'-bb')+i(ab'+a'b)
         et donc  partie réelle de x: aa'-bb'     
                  partie imaginaire de x:  ab'+a'b

A ce stade, nous obtenons:
sin 45° = 1 / 1.414213562 = 0.707107 
cos 45° = 1 / 1.414213562 = 0.707107
Multiplions maintenant le résultat par:
1+Ki avec K=2-1 (correspondant à un angle de 26.56505°). Encore une fois, 1+Ki et pas 1-Ki pour nous rapprocher des 53°.
Nous obtenons:
un angle de 45°+26.56505°=71.56505° (18.56505° en trop)
une partie réelle=1.581138830 cos 71.56505°  mais aussi = 1-K*1 = 0.5
une partie imaginaire=1.581138830 sin 71.56505° mais aussi = K*1+1 = 1.5

A ce stade, nous obtenons:
cos 71.56505° = 0.5 / 1.581138830 = 0.316228
sin 71.56505° = 1.5 / 1.581138830 = 0.948683
Multiplions maintenant le résultat par:
1-Ki avec K=2-2 (correspondant à un angle de 14.03624°). Cette fois, 1-Ki pour compenser les 18.56505° excédentaires.
Nous obtenons:
un angle de 71.56505°-14.03624°=57.52881° (4.52881° en trop)
une partie réelle=1.629800601 cos 57.52881° mais aussi = 1/2+K*3/2 = 0.875 
une partie imaginaire=1.629800601 sin 57.52881° mais aussi = -K*1/2+3/2 = 1.375

A ce stade, nous obtenons:
cos 57.52881° = 0.875 / 1.629800601 = 0.536875
sin 57.52881° = 1.375 / 1.629800601 = 0.843661
Multiplions maintenant le résultat par:
1-Ki avec K=2-3 (correspondant à un angle de 7.12502°). Cette fois, 1-Ki pour compenser les 4.52881° excédentaires.
Nous obtenons:
un angle de 57.52881°-7.12502°=50.40379° (manquent 2.59621°)
une partie réelle=1.642484066 cos 50.40379° mais aussi = 0.875+K*1.375 = 1.046875 
une partie imaginaire=1.642484066 sin 50.40379° mais aussi = -K*0.875+1.375 = 1.265625

A ce stade, nous obtenons:
cos 50.403791° = 1.046875 / 1.642484066 = 0.637373
sin 50.403791° = 1.265625 / 1.642484066 = 0.770555
Multiplions maintenant le résultat par:
1+Ki avec K=2-4 (correspondant à un angle de 3.57633°). Cette fois, 1+Ki pour compenser le manque de 2.59621°.
Nous obtenons:
un angle de 50.40379°+3.57633°=53.98012° (0.98012° en trop)
une partie réelle=1.645688916 cos 53.98012° mais aussi = 1.046875-K*1.265625 = 0.967773
une partie imaginaire=1.645688916 sin 53.98012° mais aussi = K*1.046875+1.265625 = 1.331055

A ce stade, nous obtenons:
cos 53.980126° = 1.046875 / 1.645688916 = 0.588066
sin 53.980126° = 1.265625 / 1.645688916 = 0.808813
Multiplions maintenant le résultat par:
1-Ki avec K=2-5 (correspondant à un angle de 1.78991°). Cette fois, 1-Ki pour compenser l'excédent de 0.98012°.
Nous obtenons:
un angle de 53.980126°-1.78991°=52.190215° (0.809785° trop peu)
une partie réelle=1.646492279 cos 52.190215° mais aussi = 0.967773-K*1.331055 = 1.009369
une partie imaginaire=1.646492279 sin 52.190215° mais aussi = K*0.967773+1.331055 = 1.300812

A ce stade, nous obtenons:
cos 52.190215° = 1.009369 / 1.646492279 = 0.613042
sin 52.190215° = 1.300812 / 1.646492279 = 0.790050
Multiplions maintenant le résultat par:
1+Ki avec K=2-6 (correspondant à un angle de 0.89517°). Cette fois, 1+Ki pour compenser le manque de 0.809785°.
Nous obtenons:
un angle de 52.190215°+0.89517°=53.085389° (0.085389° en trop)
une partie réelle=1.646693254  cos 53.085389° mais aussi = 1.009369-K*1.300812 = 0.989044
une partie imaginaire=1.646693254 sin 53.085389° mais aussi = K*1.009369+1.300812 = 1.316583

A ce stade, nous obtenons:
cos 53.085389° = 0.989044 / 1.646693254 = 0.600624
sin 53.085389° = 1.316583 / 1.646693254 = 0.799532
Si nous arrêtons ici, nous obtenons:
cos 53.085389° = 0.600624 alors que cos 53° = 0.60181502315204827991797700044149
sin 53.085389° = 0.799532 alors que sin 53° = 0.79863551004729284628400080406894
Si nous avions continué, nous aurions obtenu:
cos 52.999719° = 0.601819 et donc une erreur <5*10-6
sin 52.999719° = 0.798633 et donc une erreur <5*10-6

Il est évident que l'explication qui précède est 'brute'. En effet, pour calculer sin 3°, il serait stupide de passer par 
45°, 45°-26.56505°=18.43495°,18.43495°-14.03624°=4.39871°, 4.39871°-7.12502°=-2.72631°
alors qu'il est possible de commencer par 3.576330° qui nous met à 0.5° de l'angle cherché.

De même pour un angle supérieur à 90°, on ne commence pas avec 0° mais avec 90° (et donc 0+i 1).

Voici une implémentation en Turbo Pascal (dans l'archive, vous retrouvez la version 'brute' et celle ci-dessous):
PROGRAM Cordic2;

USES CRT;

VAR Angle,a,PartR,PartI,K,Tmp,Gain,Diff:DOUBLE;
    Angles,Gains,Puiss2:ARRAY[0..50] OF DOUBLE;
    i,j:BYTE;
    Signe:SHORTINT;

PROCEDURE Init;
VAR i:BYTE;
    K:REAL;
BEGIN
K:=2;
FOR i:=0 TO 50 DO
  BEGIN
  K:=K/2;
  Puiss2[i]:=K;
  Angles[i]:=ARCTAN(K)*180/Pi;
  Gains[i]:=SQRT(1+K*K);
  END;
END;

PROCEDURE ReduitAngle(VAR a:DOUBLE);
BEGIN
IF a>0 THEN WHILE NOT (a<360) DO
              a:=a-360
ELSE WHILE NOT (a>0) DO
       a:=a+360
END;

FUNCTION Search(D:DOUBLE;i:BYTE):BYTE;
VAR j:BYTE;
BEGIN
j:=i;
WHILE NOT (ABS(D-Angles[j])<D) DO
  INC(j);
Search:=j
END;

BEGIN
Init;
WRITE('Introduisez l''angle : ');READLN(a);
ReduitAngle(a);
i:=0;Gain:=1;
IF a>90 THEN BEGIN
             PartR:=0;PartI:=1;Angle:=90;
             END
ELSE BEGIN
     PartR:=1;PartI:=0;Angle:=0;
     END;
Diff:=a-Angle;
WHILE NOT (ABS(Diff)<1E-3) DO
  BEGIN
  j:=Search(ABS(Diff),i);
  IF Diff<0 THEN Signe:=-1
  ELSE Signe:=+1;
  Angle:=Angle+Signe*Angles[j];
  Diff:=a-Angle;
  K:=Puiss2[j];
  Gain:=Gain*Gains[j];
  Tmp:=PartR-Signe*K*PartI;
  PartI:=Signe*K*PartR+PartI;
  PartR:=Tmp;
  WRITELN(Angle:8:6,' sinus=',PartI/Gain:8:6,' & cosinus=',PartR/Gain:8:6);
  i:=j;
  READKEY;
  END;
WRITELN;
WRITELN('sinus=',PartI/Gain:8:6,' & cosinus=',PartR/Gain:8:6);
END.
Il est aussi certain qu'il est possible de se passer des tableaux et de ne calculer les angles et gains que quand ils sont utiles.

Page précédente.
Page d'accueil.