# Very High Speed Integrated Circuit Hardware Description Language/Réalisation d'un coprocesseur CORDIC

Nous appellerons un coprocesseur dans la suite de ce cours tout périphérique suffisamment sophistiqué pour :
* faire des calculs arithmétiques au sens large
* nécessiter un processeur pour être testé
* éventuellement être capable d'exécuter un programme
Les GPU (des cartes graphiques) se rangent dans cette catégorie mais aussi les processeurs mathématiques (Unités de calcul en virgule flottante (ou FPU)) capables de calculer des fonctions simples ou compliquées comme les fonctions trigonométriques ou les fonctions hyperboliques.

## Introduction à CORDIC

L'algorithme CORDIC est un algorithme itératif permettant de faire des calculs de fonctions trigonométriques, de fonctions hyperboliques ou de fonctions linéaires. Cet algorithme est difficile à appréhender lors d'une première approche. Heureusement, utiliser un algorithme ne nécessite pas de comprendre comment il fonctionne dans les détails. Mais sa modification ou son adaptation par contre le nécessite.

### Domaines d'utilisation de l'algorithme CORDIC

L'algorithme CORDIC peut naturellement être simplement utilisé pour des calculs. Pourtant, dans un livre sur les processeurs enfouis, nous ne resterons pas dans le domaine de l'abstraction. Nous allons donc lui consacrer au moins deux domaines d’applications :
* robotique mobile : calcul des coordonnées x et y d'un robot se déplaçant à l'aide de deux moteurs à courant continu. Lisez AVR et robotique : ASURO dans un autre projet et particulièrement la section "Marier l'électronique et la cinématique" vous y trouverez des formules faisant intervenir des sinus et cosinus.
* commande de moteur synchrone

### Propriétés des fonctions trigonométriques

Les calculs se font avec des nombres. Puisqu'ici nous désirons étudier des "sinus" et "cosinus" et que ces nombres sont entre 0 et 1 il nous faudra réfléchir sur la façon de représenter ces nombres. C'est l’objet de la prochaine section.

## Représentation des nombres à virgule

Il existe au moins deux méthodes pour faire cela :
* la représentation virgule flottante
* la représentation virgule fixe
Nous allons présenter les deux méthodes même si nous avons l'intention de nous attarder beaucoup plus sur la virgule fixe.

### La représentation virgule flottante

L'article Virgule_flottante de wikipédia est suffisamment complet pour ne pas être repris ici. Il peut d'ailleurs être complété par la lecture de l’article détaillé sur la norme IEEE 754 correspondante.
Voici en figure récapitulative le codage des nombres flottants.
Il existe essentiellement deux formats :
| | Encodage | Signe | Exposant | Mantisse | Valeur d'un nombre | Précision | Chiffres significatifs |
| Simple précision | 32 bits | 1 bit | 8 bits | 23 bits | $(-1)^{S}\times M\times 2^{(E-127)}$ | 24 bits | environ 7 |
| Double précision | 64 bits | 1 bit | 11 bits | 52 bits | $(-1)^{S}\times M\times 2^{(E-1023)}$ | 53 bits | environ 16 |
Sur ces deux formats, seul le format simple précision sur 32 bits nous intéresse. Il peut se résumer avec la figure :
Pour un nombre normalisé, le bit de poids fort de la mantisse est toujours à 1 et il n'est alors pas représenté. Il est alors qualifié de bit implicite. La conséquence est que l'expression du tableau pour la représentation 32 bits doit être comprise comme :
* une différence entre le M de la formule et la Mantisse : M = 1,Mantisse
* la valeur du nombre est donnée par $(-1)^{S}\times 1,Mantisse\times 2^{(E-127)}$
Si le bit de poids fort de la mantisse (ici bit 23 implicite) est nul (en général pour un exposant décalé nul) alors le nombre est dit dénormalisé.

### La représentation virgule fixe

Le format Virgule fixe n’est pas normalisé contrairement au format virgule flottante. Le seul point commun que l’on retrouve dans tous les formats virgule fixe est qu’il est signé. Il peut donc représenter des nombres positifs ou négatifs en complément à 2. La virgule étant placée par défaut à un endroit fixe on utilise en général une notation du genre Qm.n pour le désigner où m représente le nombre de bits avant la virgule et n le nombre de bit après la virgule et donc m+n est la taille de la représentation.
Par exemple nous allons utiliser un format Q3.13 dans la suite de ce chapitre.
| $b_{15}$ | $b_{14}$ | $b_{13}$ | $b_{12}$ | $b_{11}$ | $b_{10}$ | $b_{9}$ | $b_{8}$ | $b_{7}$ | $b_{6}$ | $b_{5}$ | $b_{4}$ | $b_{3}$ | $b_{2}$ | $b_{1}$ | $b_{0}$ |
| S | $2^{1}$ | $2^{0}$ | $2^{-1}$ | $2^{-2}$ | $2^{-3}$ | $2^{-4}$ | $2^{-5}$ | $2^{-6}$ | $2^{-7}$ | $2^{-8}$ | $2^{-9}$ | $2^{-10}$ | $2^{-11}$ | $2^{-12}$ | $2^{-13}$ |
Il s'agit d'un format 16 bits (3+13) et adapté à la trigonométrie en radian.
Pour manipuler ce format nous proposons quelques utilitaires en C. Ils sont basés sur le fait que la librairie standard du compilateur C sait manipuler les nombres à virgule à partir du moment où ils sont en virgule flottante. La lecture à partir du clavier ainsi que l’affichage sur un écran de ces nombres est aussi dans la librairie standard.
Nous proposons donc deux routines permettant les transformations de notre format virgule fixe Q3.13 vers le format flottant.

#### Format Q3.13 vers flottant

Nous vous proposons deux sous-programmes destinés à une conversion du format Q3.13 vers le format float du C. Nous avons écrit ces sous-programme avec l'intention un jour de les implanter en matériel. C'est pour cela qu'ils peuvent vous paraître obscurs. En effet une simple division flottante peut naturellement faire l'affaire.
Sans entrer dans les détails, nous vous proposons des sous-programmes réservés aux architectures de type PC :
* Une simple division suffit à faire cette conversion : il suffit de diviser par 1 (soit 0x2000 en Q3.13). Nous avons choisi ces méthodes un peu complexe avec l'idée derrière la tête, de les matérialiser un jour
La deuxième routine n'utilisant pas "pow" n'a pas beaucoup d'intérêt pour un PC. Mais elle a été développée en vue d'un portage sur AVR, ce qui a été fait quelques années après (2019).
Comme nous l'avons énoncé au tout début de cette section, une simple routine universelle pourrait donc être réalisée par :
Attachons-nous maintenant à la transformation inverse.

#### Format flottant vers Q3.13

La transformation du format flottant vers le Q3.13 se fait par :
Cette première routine ne contient pas de vérification d'erreur tandis que la deuxième affiche une erreur mais retourne quand même un nombre. Le format "float" du c étant plus précis que notre format Q3.13, cette routine peut échouer. Elle reste cependant utile dans le cadre de CORDIC.

#### Tests des deux routines précédentes avec un Arduino

Pour tester le bon fonctionnement des transformations de formats, rien de tel qu'un programme qui transforme dans les deux sens pour constater si oui ou non on retrouve le même résultat. Voici donc un exemple sous Arduino.

## Un peu de théorie sur CORDIC

Dans cette section, nous allons détailler l'algorithme CORDIC. Nous présentons un schéma simplifié de principe de calcul où vous voyez apparaître trois entrées $(x_{0},y_{0},z_{0})$ et trois sorties $(x_{n},y_{n},z_{n})$. Ce schéma est simplifié dans le sens où l’on a omis l'horloge en entrée (il faut bien une horloge pour faire avancer un calcul).
Les données à gauche représentent l'initialisation de l'algorithme. La donnée essentielle est $z_{0}$ qui représente un angle qui pour nous sera exprimé en radians. Les autres données seront prises comme $x_{0}=1$ et $y_{0}=0$ dans un premier temps pour expliquer l'algorithme.
Les sorties à droites représenteront le sinus et le cosinus ainsi que l'erreur réalisée sur l'angle.
Qu'y a-t-il dans le bloc ? Un ensemble de trois équations de récurrences que nous allons examiner maintenant.

### Des équations de récurrences simples

Nous avons au départ trois équations de récurrences simples :
$\left\{{\begin{matrix}x_{i+1}=x_{i}-\sigma _{i}y_{i}\delta _{i}\\y_{i+1}=y_{i}+\sigma _{i}x_{i}\delta _{i}\\z_{i+1}=z_{i}-\sigma _{i}\alpha _{i}\end{matrix}}\right.$
avec $\sigma _{i}\in \{-1;1\}$.
On expliquera plus tard les relations
$\delta _{i}=2^{-i}$ et $\alpha _{i}=atan(2^{-i})$ (constantes prédéfinies).
Ce que font ces trois équations est représenté sur la figure ci-contre où les variables $z_{i}$ ont été remplacées par une notation plus conventionnelle $\alpha _{i}$ pour des angles. Nous allons expliciter cette figure maintenant.
Dans cette figure vous voyez apparaître un ensemble de rotations qui aboutissent à l'angle dont on cherche le sinus et le cosinus (vecteur rouge de la figure). Vous avez certainement remarqué que les rotations ne sont pas des vraies rotations (qui conservent les longueur) puisqu'elles sont systématiquement associées à une dilatation. La dilatation dépend de l'angle mais :
* Si l’on appelle R, K et H (pour hypoténuse) les côtés du triangle, on a les relations :
** $K=R\cdot tan(\alpha )$ donc $H^{2}=R^{2}+R^{2}\cdot tan(\alpha )$ et ainsi
** $H={{R} \over {cos(\alpha )}}$ (relation importante qui sera retrouvée autrement plus tard).
Ce qui n’apparaît pas sur la figure est que les angles sont déterminés à l'avance (on reviendra sur cette propriété plus tard).
* Quelle que soit la séquence de rotations exécutées, si elles sont assez nombreuses la dilatation finale sera toujours la même.
Ces trois propriétés justifient la méthode CORDIC. Mais regardons d'un peu de plus près les équations de récurrence et donnons-en quelques caractéristiques à travers des remarques.
$\sigma _{i}y_{i}\delta _{i}$
La multiplication par $\sigma _{i}$ n’est pas une multiplication car cette valeur est soit +1 (rien à faire) soit -1 (un complément à deux à réaliser). Or le complément à deux se réalise avec des inverseurs et une addition de 1 ! donc pas de multiplication.
$\delta _{i}$
Maintenant que tous les lecteurs ont entraperçu une relation entre CORDIC et les rotations, nous allons examiner cela de manière un peu plus formelle.

### Relation avec les rotations

Cette partie s'inspire un peu d'un sugjet d'agrégation Génie électrique (2006) : sujets agrégation. Les sujets d'agrégation ne sont pas soumis au copyright. De toute façon ce sujet d'agrégation a été fortement remanié pour le rendre plus accessible. Il en résulte un ensemble de questions avec ses solutions.

#### Écriture matricielle de CORDIC

Écrire les deux premières équations sous la forme matricielle en calculant la matrice Ci : $v_{i+1}=C_{i}\cdot v_{i}$
avec $v_{i}=(x_{i},y_{i})^{T}$ et $v_{i+1}=(x_{i+1},y_{i+1})^{T}$. (On note donc tout simplement l'entrée et la sortie sous forme de deux composantes vectorielles pour pouvoir introduire cette matrice).
$C_{i}={\begin{pmatrix}1&-\sigma _{i}\delta _{i}\\[3pt]\sigma _{i}\delta _{i}&1\\\end{pmatrix}}$

#### Comment introduire les matrices de rotation

La figure montre clairement un ensemble de rotations. Rappelons que les rotations dans un plan sont décrites par des matrices de rotation rappelées ci-après : $R(\alpha )={\begin{pmatrix}\cos \alpha &-\sin \alpha \\[3pt]\sin \alpha &\cos \alpha \\\end{pmatrix}}$ (rotation d'angle α)
Cette matrice fait tourner le plan d'un angle α.
Elle peut être généralisée sous la forme indicée :
$R(\alpha _{i})={\begin{pmatrix}\cos \alpha _{i}&-\sigma _{i}\sin \alpha _{i}\\[3pt]\sigma _{i}\sin \alpha _{i}&\cos \alpha _{i}\\\end{pmatrix}}$
où le paramètre $\sigma _{i}$ pour décrire le sens de rotation ($\sigma _{i}=+1$ rotation dans le sens trigonométrique).
L'idée qui vient à l'esprit est donc de tenter d'écrire la matrice $C_{i}$ de la section précédente à l'aide d'une matrice de rotation
Écrire la matrice $C_{i}$ sous la forme $C_{i}=K_{i}(\alpha _{i})\cdot R_{i}(\alpha _{i})$ en prenant soin d'expliciter $K_{i}(\alpha _{i})$ et finalement $\delta _{i}$.
$C_{i}={\begin{pmatrix}1&-\sigma _{i}\delta _{i}\\[3pt]\sigma _{i}\delta _{i}&1\\\end{pmatrix}}={1 \over cos(\alpha _{i})}\cdot {\begin{pmatrix}\cos \alpha _{i}&-\sigma _{i}\sin \alpha _{i}\\[3pt]\sigma _{i}\sin \alpha _{i}&\cos \alpha _{i}\\\end{pmatrix}}$
qui donne immédiatement $k_{i}(\alpha _{i})={1 \over cos(\alpha _{i})}$ et $\delta _{i}$ devient : $\delta _{i}={sin(\alpha _{i}) \over cos(\alpha _{i})}$
Prenons du recul par rapport aux équations : en clair les relations de récurrences CORDIC réalisées par la matrice $C_{i}$ sont équivalentes à une rotation $R(\alpha _{i})$ multipliée par $1 \over cos(\alpha _{i})$ qui est plus grand que 1. C'est donc équivalent à une rotation et une dilatation.
Tiens tiens, n'a-t-on pas déjà évoqué cette propriété dans ce chapitre d'une dilatation par ${1 \over cos(\alpha _{i})}$ ?
Arrivé à ce point, nous espérons que le lecteur a compris (ou au moins entrevu) la relation entre les équations de récurrences CORDIC et les rotations. Rappelons que l’on cherche à calculer un cosinus et un sinus. Il reste cependant des détails à régler pour faire fonctionner CORDIC.

#### Calculer les équations de récurrences

Soit $v_{0}$ la valeur initiale de $v_{i}$, exprimer $v_{i}$ en fonction de $v_{0}$, $C_{i-1}$, $C_{i-2}$, ... puis en fonction de $v_{0}$, $R_{i-1}$, $R_{i-2}$, $K_{i-1}$, $K_{i-1}$ ...
$v_{i}=C_{i-1}\cdot C_{i-2}\cdot ...\cdot C_{0}\cdot v_{0}$
$v_{i}=K_{i-1}\cdot R_{i-1}\cdot K_{i-2}\cdot R_{i-2}\cdot ...\cdot K_{0}\cdot R_{0}\cdot v_{0}$
Il est temps de travailler sur la troisième équation de récurrence que nous avons laissé tombé jusqu'à présent.

#### Travail sur les angles

Premier détail : nous allons prendre à partir de maintenant $\delta _{i}=2^{-i}$. C'est cette subtilité qui fait toute la valeur de l'algorithme CORDIC. Pourquoi ? Parce que ce terme $\delta _{i}$ apparaît toujours comme terme multiplicatif dans les équations de récurrences :
$\left\{{\begin{matrix}x_{i+1}=x_{i}-\sigma _{i}y_{i}\delta _{i}\\y_{i+1}=y_{i}+\sigma _{i}x_{i}\delta _{i}\\z_{i+1}=z_{i}-\sigma _{i}\alpha _{i}\end{matrix}}\right.$
et qu'une propriété intéressante d'une multiplication par une puissance de 2 est qu'elle peut être remplacée par un simple décalage. Autrement dit choisir $\delta _{i}=2^{-i}$ va permettre d'économiser des multiplications (qui est un opérateur arithmétique considéré comme complexe, lent et cher).
$\alpha _{i}=atan(2^{-i})$
Calculer $\alpha _{i}$ pour i=0, 1, 2, ..., 13 si, comme on le rappelle, $\delta _{i}={sin(\alpha _{i}) \over cos(\alpha _{i})}$ (question précédente).
Réponse : On donne les valeurs précalculées :
| i | α$_i$ = tan$^{−1}$ (2$^{−i}$) | α$_i$ = tan$^{−1}$ (2$^{−i}$) |
| i | Degrees | Radians |
| 0 | 45.00 | 0.7854 |
| 1 | 26.57 | 0.4636 |
| 2 | 14.04 | 0.2450 |
| 3 | 7.13 | 0.1244 |
| 4 | 3.58 | 0.0624 |
| 5 | 1.79 | 0.0312 |
| 6 | 0.90 | 0.0160 |
| 7 | 0.45 | 0.0080 |
| 8 | 0.2238 | 0.0039 |
| 9 | 0.1119 | 0.0019 |
| 10 | 0.05595 | 0.00098 |
| 11 | 0.02798 | 0.00049 |
| 12 | 0.01399 | 0.00024 |
| 13 | 0.00699 | 0.00012 |
Deuxième détail : nous avons laissé de côté $\sigma _{i}$ qui est un signe + ou -. Il est grand temps de voir comment il est calculé. C'est tout simplement le signe de l'angle calculé : $\sigma _{i}=signe(z_{i})$. Ce signe $\sigma _{i}$ est connu (puisque $z_{i}$ l'est) lors du calcul. Le calcul précédent peut être réalisé tel quel dans du matériel (qui peut calculer en parallèle).

#### Le problème du gain

Depuis le début de ce chapitre nous savons qu'un problème de gain était à prendre en compte. En effet, nous avons montré très tôt que l'algorithme réalisait en fait une rotation $R_{z}(\alpha _{i})$ mais multiplié par $K_{i}={1 \over {cos(\alpha _{i})}}$. Pour avoir un cosinus et un sinus il nous faudra donc calculer $K=\prod _{i=0}^{n-1}K_{i}={\prod _{i=0}^{n-1}{1 \over {cos(\alpha _{i})}}}$ et multiplier le résultat final par cette valeur.
On peut éviter cependant cette multiplication si au lieu d'initialiser $x_{0}$ à 1, on réalise : $x_{0}={1 \over {\prod _{i=0}^{n-1}{1 \over {cos(\alpha _{i})}}}}={1 \over K}$.
Les valeurs numériques correspondantes sont facilement évaluées avec un tableur ($K=\prod _{i=0}^{n-1}K_{i}=1,64676$ d'inverse 0,60725).
Voila. Maintenant tous les détails de calculs sont en place. Pour mettre cela en œuvre dans un premier temps, nous allons utiliser un tableur et regarder CORDIC fonctionner puis nous ferons la même chose en langage C pur.

## Utiliser un tableur pour faire fonctionner CORDIC

Nous nous proposons dans cette section d’utiliser le tableur LibreOffice pour examiner comment fonctionne CORDIC. Nous pensons, même si nous n'avons pas essayé, qu'EXEL utiliserait les mêmes primitives.
Nous allons décomposer le travail en deux parties.

### Faire fonctionner la troisième équation de récurrence

La troisième équation de récurrence est un calcul d'angle $z_{i+1}$ en fonction de $z_{i}$. Comme nous l'avons déjà affirmé, cette équation peut fonctionner seule. Voici comment on peut procéder avec un tableur (pour lequel les colonnes sont des lettres et les lignes des chiffres, format par défaut de LibreOffice):
| | A | B | C | D | E | F |
| 1 | i | delta_i | alpha_i | K_i | z_i | sigma_i |
| 2 | 0 | =PUISSANCE(2;-A2) | =ATAN(B2) | =1/COS(C2) | =PI/3 | =SI(E2<0;-1;1) |
| 3 | 1 | =PUISSANCE(2;-A3) | =ATAN(B3) | =1/COS(C3) | =E2-C2*F2 | =SI(E3<0;-1;1) |
que l’on peut étendre vers le bas.
Le point essentiel du tableau est la case E2 (qui contient =PI/3) qui est l'angle duquel on cherche le sinus et le cosinus. C'est la seule case à changer pour faire le calcul sur un autre angle.
L'intérêt d’utiliser un tableur est de voir comment varie cette variable z_i avec les itérations : elle converge toujours vers 0 !
Voici le résultat complet :
Ce graphe montre une convergence vers 0 qu’il n’est pas inutile de retenir pour comprendre les différentes utilisations de CORDIC évoquées plus loin.

### Faire fonctionner l'ensemble

On complète le tableau par les deux colonnes représentant les deux équations manquantes (colonnes G et H) :
| | A | B | C | D | E | F | G | H | I |
| 1 | i | delta_i | alpha_i | K_i | z_i | sigma_i | x_i | y_i | angle_i |
| 2 | 0 | =PUISSANCE(2;-A2) | =ATAN(B2) | =1/COS(C2) | =PI/3 | =SI(E2<0;-1;1) | 0,6073 | 0 | 0 |
| 3 | 1 | =PUISSANCE(2;-A3) | =ATAN(B3) | =1/COS(C3) | =E2-C2*F2 | =SI(E3<0;-1;1) | =G2-F2*H2*B2 | = H2+G2*F2*B2 | =I2+F2*C2 |
On a ajouté pour information une colonne à droite qui calcule la somme des rotations pour montrer qu'elle converge bien vers l'angle souhaité.
Voici le résultat complet :
Encore une fois, seule la case E2 est à changer pour réaliser un autre calcul.
$-\pi \over 2$
$+\pi \over 2$
* La troisième équation de récurrence peut être indistinctement celle sur $z_{i}$ ou celle sur $angle_{i}$. Mais alors le calcul du signe $sigma_{i}$ se fait différemment. Les dessins d'introduction montraient plutôt un angle qui convergeait vers l'angle final alors que notre premier algorithme CORDIC dans le tableur utilise $z_{i}$.

## Utiliser et adapter l'algorithme de Wikipédia

L'article CORDIC dans wikipédia propose un algorithme en C (enfin proposait : voir remarque ci-dessous). Prenez-le avec un copier coller et à l'aide d'un éditeur faites-en un fichier que l’on peut appeler "cordicWiki.c" par exemple. Sa compilation sous Linux ainsi que son exécution peut se faire par :
$-\pi \over 2$
$+\pi \over 2$
L'utilisation en cours de wikipédia nous amène à faire la remarque suivante :
Pour pallier à cet inconvénient, nous le remettons ici :
Le fonctionnement de cet algorithme semble parfait. Pourtant, nous allons essayer de l'améliorer pour un fonctionnement dans une architecture bien plus petite qu'un PC. En effet, compilé tel quel avec le compilateur avr-gcc, il faut plus de 8 ko de mémoire, ce qui n’est pas toujours possible pour certains processeurs (ATMega8 et ATTiny861) que nous utilisons dans ce cours.

### Analyse des ressources calculs nécessaires

Une lecture de ce programme montre qu’il utilise le type "double" qui est le type flottant sur 64 bits. C'est un type peu adapté aux architectures 8 bits que nous désirons cibler.
Ce programme fonctionne avec les librairies de calcul en double (librairie mathématique qu'utilisent pow(2,-i) et atan(Pow2)) et nous désirons à tout prix l'éviter.

### Simplifier le programme

La technique habituelle pour faire cette simplification est de choisir un format virgule fixe en lieu et place des variables de type "double". Nous allons utiliser le format Q3.13 qui veut dire que trois bits de poids forts sont la partie entière et les treize bits restants sont la partie fractionnaire.
D'autre part les multiplications par pow(2,-i) (soit avec des notations plus habituelles $\delta _{i}=2^{-i}$) sont remplacées par de simples décalages comme cela a déjà été évoqué.
Enfin le calcul atan(Pow2) (soit avec des notations plus habituelles $arctan(\delta _{i})=arctan(2^{-i})$) sera remplacé par un tableau de valeurs pré-calculées.
Ce travail a déjà été fait dans ICI (dans un autre livre) où vous avez un programme complet sans aucune multiplication.

### Problème d'arrondi

Sans entrer dans les détails, notons que cet algorithme peut être amélioré. Toute représentation des nombres présente des imperfections. Ici, supposons que l’on ait 0xFFFF en format Q3.13. Un décalage vers la droite le transformerait en 0xFFFF alors qu’il faudrait le transformer en 0x0000 !