Très probablement, vous connaissez les ratios suivants de l'école:

$$

$$\ sin (\ alpha + \ beta) = \ sin \ alpha \ times \ cos \ beta + \ cos \ alpha \ times \ sin \ beta \\ \ cos (\ alpha + \ beta) = \ cos \ alpha \ times \ cos \ beta - \ sin \ alpha \ times \ sin \ beta$$

Lorsque vous vous êtes familiarisé avec cette formule pour la première fois dans votre enfance, il est fort probable que votre premier sentiment ait été la douleur, car il faut se souvenir de cette formule. C'est très mauvais, car en fait vous n'avez pas besoin de vous souvenir de cette formule - elle s'affiche lorsque vous faites pivoter le triangle sur du papier. En fait, je fais de même lorsque j'écris cette formule. Cette interprétation apparaîtra au milieu de cet article. Mais maintenant, pour laisser tout le plaisir pour plus tard et reporter le moment où vous dites «Eureka!», Réfléchissons à la raison pour laquelle nous devrions même penser à cette formule.

Entrée

Les fonctions trigonométriques sin() et cos() probablement les plus populaires en infographie, car elles sont la base pour décrire n'importe quelle forme ronde de manière paramétrique. Parmi les lieux de leur application possible figurent la génération de cercles ou d'objets volumétriques de révolution, lors du calcul de la transformée de Fourier, la génération procédurale d'ondes sur le plan de l'eau, des générateurs pour un logiciel de synthèse sonore, etc. Dans tous ces cas, sin() et cos() sont appelés à l'intérieur de la boucle, comme ici:

 for(int n=0; n < num; n++) { const float t = 2.0f*PI*(float)n/(float)num; const float s = sinf(t); const float c = cosf(t); // do something with s and c ... } 

Nous commençons à réécrire le cycle de manière incrémentielle (voir le code ci-dessous), il est donc plus facile pour nous d'imaginer qu'à l'itération n ce cycle avec la phase t , l'itération suivante, n+1 , considérera sin() et cos() pour t+f . En d'autres termes, sin(t) et cos(t) sont comptés pour nous et nous devons compter sin(t+f) et cos(t+f) :

 const float f = 2.0f*PI/(float)num; const float t = 0.0f; for( int n=0; n < num; n++ ) { const float s = sinf(t); const float c = cosf(t); // do something with s and c ... t += f; } 

Peu importe comment nous avons obtenu t et quelle est sa plage de valeurs (dans l'exemple ci-dessus - $$$[0; 2 \ pi]$) La seule chose qui nous dérange, c'est qu'il existe une boucle qui appelle constamment sin() et cos() avec un paramètre qui augmente par pas constants (dans ce cas, $$$\ frac {2 \ pi} {\ text {num}}$) Cet article explique comment optimiser ce code pour la vitesse de telle sorte que les mêmes calculs puissent être effectués sans utiliser les fonctions sin() ou cos() (dans la boucle interne), ni même la fonction plus rapide sincos() .

Mais si vous regardez la première formule de l'article, nous verrons que si $$$f = \ alpha$et $$$t = \ beta$nous pouvons réécrire ceci

 sin(t+f) = sin(f)*cos(t) + cos(f)*sin(t) cos(t+f) = cos(f)*cos(t) - sin(f)*sin(t) 

ou en d'autres termes

 new_s = sin(f)*old_c + cos(f)*old_s new_c = cos(f)*old_c - sin(f)*old_s 

Puisque f est une constante, sin(f) et cos(f) également. Appelez-les respectivement a et b :

 new_s = b*old_c + a*old_s new_c = a*old_c - b*old_s 

Cette expression peut être directement utilisée dans notre code, puis nous obtenons une boucle dans laquelle des fonctions trigonométriques coûteuses (et même pas) sont appelées!

 const float f = 2.0f*PI/(float)num; const float a = cosf(f); const float b = sinf(f); float s = 0.0f; float c = 1.0f; for( int n=0; n < num; n++ ) { // do something with s and c ... const float ns = b*c + a*s; const float nc = a*c - b*s; c = nc; s = ns; } 

Interprétation

Jusqu'à présent, nous avons joué aveuglément avec les mathématiques, sans comprendre ce qui se passe réellement. Réécrivons la boucle intérieure comme ceci:

$$

$$s_ {n + 1} = s_na + c_nb \\ c_ {n + 1} = c_na - s_nb$$

Certains d'entre vous ont peut-être remarqué qu'il s'agit d'une formule pour faire tourner un objet dans un espace à deux dimensions. Si vous ne comprenez toujours pas cela, peut-être que la forme matricielle vous aidera.

$$

$$\ left (\ begin {matrix} s_ {n + 1} \\ c_ {n + 1} \ end {matrix} \ right) = \ left (\ begin {matrix} a & b \\ -b & a \ end {matrix} \ right) \ cdot \ left (\ begin {matrix} s_ {n} \\ c_ {n} \ end {matrix} \ right)$$

En fait, sin(t) et cos(t) peuvent être regroupés en un vecteur de longueur 1 et dessinés comme un point sur l'écran. Appelez ce vecteur x . Alors $$$x = \ {\ sin \ beta, \ cos \ beta \}$ . Ainsi, la forme d'expression vectorielle est $$$x_ {n + 1} = Rx_n$où $$$R = \ left (\ begin {matrix} a & b \\ - b & a \ end {matrix} \ right)$ . On voit que notre boucle effectue une petite rotation bidimensionnelle à chaque itération pour que x tourne en cercle lors de l'exécution du cycle. Chaque itération x tourne de $$$\ frac {2 \ pi} {\ text {num}}$radian.
Donc en gros

$$

$$\ sin (\ alpha + \ beta) = \ sin \ alpha \ times \ cos \ beta + \ cos \ alpha \ times \ sin \ beta \\ \ cos (\ alpha + \ beta) = \ cos \ alpha \ times \ cos \ beta - \ sin \ alpha \ times \ sin \ beta$$

c'est la formule du mouvement ponctuel $$$x = \ {\ cos \ alpha, \ sin \ alpha \}$ circonférence par incréments de $$$\ beta$radian. Pour ce faire, nous allons construire l'un des deux axes de rotation, par exemple, $$$u = \ {\ cos \ beta, \ sin \ beta \}$ . La première composante de la rotation est la projection. $$$x$sur $$$u$. Depuis $$$x$et $$$u$normalisé (avoir une longueur de 1), la projection est leur produit scalaire. Par conséquent $$$s = x \ cdot u = \ sin \ alpha \ cdot \ cos \ beta + \ cos \ alpha \ cdot \ sin \ beta$, et bien sûr le deuxième composant est l'anti-projection, qui peut être trouvée en se projetant sur l'axe perpendiculaire, $$$v$. Nous pouvons créer ce vecteur en développant les coordonnées $$$u$et changez le signe à l'opposé à la première coordonnée: $$$c = x \ cdot v = \ cos \ alpha \ cdot \ cos \ beta + \ sin \ alpha \ cdot \ sin \ beta$

Remarques

Habituellement, vous devriez pouvoir effectuer ces minuscules rotations encore et encore. En fait $$$| R | = \ left | \ begin {matrix} a & b \\ - b & a \ end {matrix} \ right | = a ^ 2 + b ^ 2 = \ sin ^ 2 \ alpha + \ cos ^ 2 \ alpha = 1$ , ce qui signifie que la matrice $$$R$n'augmente ni ne diminue l'espace auquel il est appliqué, ce qui signifie que $$$x$se déplacera dans un cercle parfait. Cependant, comme les ordinateurs ne sont pas précis, $$$x$se déplacera en spirale et coïncidera éventuellement avec le centre du cercle de rotation. Je n'ai pas eu de tels problèmes, mais je pense qu'ils peuvent survenir avec un très grand nombre, c'est-à-dire petits angles de braquage.

Exemple

Dans Kindercrasher , la démo de 4096 octets de 2008 (capture d'écran sur KDPV), un groupe de sphères vibre au rythme de la musique. Pour cela, j'ai compté la transformée de Fourier du son. Il existe des algorithmes qui le font en temps réel, par exemple, la FFT . Cependant, mon code devrait tenir dans quelques kilo-octets, et j'ai décidé d'aller dans l'autre sens. Au lieu d'implémenter FFT, j'ai écrit DFT par sa définition simple. Vous pouvez le vérifier sur wikipedia.

$$

$$X_k = \ sum_ {n = 0} ^ {N-1} {x_ne ^ {- \ frac {2 \ pi i} {N} kn}} \ quad k = 0,1, \ ldots, N-1$$

Ma fonction prend également un tampon audio stéréo 16 bits, x , et renvoie les 128 premières fréquences du spectre sonore du son y . Voyez comment la boucle interne est organisée, celle qui s'exécute 4096 fois: pas un seul appel aux fonctions sin() ou cos() , bien que dans d'autres implémentations, ces appels le seront.

 #include <math.h> void iqDFT12( float *y, const short *x ) { for( int i=0; i<128; i++ ) { const float wi = (float)i*(2.0f*3.1415927f/4096.0f); const float sii = sinf( wi ); const float coi = cosf( wi ); float co = 1.0f; float si = 0.0f; float acco = 0.0f; float acsi = 0.0f; for( int j=0; j<4096; j++ ) { const float f = (float)(x[2*j+0]+x[2*j+1]); const float oco = co; acco += co*f; co = co*coi - si*sii; acsi += si*f; si = si*coi + oco*sii; } y[i] = sqrtf(acco*acco+acsi*acsi)*(1.0f/32767.0f); } }