Provavelmente, você conhece as seguintes proporções da escola:

$$

$$\ 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$$

Quando você se familiarizou com esta fórmula na infância, muito provavelmente, seu primeiro sentimento foi de dor devido ao fato de que essa fórmula deve ser lembrada. Isso é muito ruim, porque na verdade você não precisa se lembrar dessa fórmula - ela é exibida quando você gira o triângulo no papel. De fato, faço o mesmo quando escrevo esta fórmula. Essa interpretação será aparente no meio deste artigo. Mas agora, para deixar toda a diversão para mais tarde e adiar o momento em que você diz “Eureka!”, Vamos pensar sobre por que deveríamos pensar sobre essa fórmula.

Entrada

As funções trigonométricas sin() e cos() provavelmente as mais populares em computação gráfica, pois são a base para descrever qualquer forma redonda de maneira paramétrica. Entre os locais de sua possível aplicação estão a geração de círculos ou objetos volumétricos de revolução, ao calcular a transformada de Fourier, a geração procedural de ondas no plano aquático, geradores para um sintetizador de som de software e assim por diante. Em todos esses casos, sin() e cos() são chamados dentro do loop, como aqui:

 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 ... } 

Começamos a reescrever o ciclo de forma incremental (veja o código abaixo), portanto é mais fácil imaginar que na iteração n deste ciclo com a fase t , a próxima iteração, n+1 , considerará sin() e cos() para t+f . Em outras palavras, sin(t) e cos(t) são contados para nós e precisamos contar sin(t+f) e 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; } 

Não importa como obtivemos t qual é seu intervalo de valores (no exemplo acima - $$$[0; 2 \ pi]$) A única coisa que nos incomoda é que existe um loop que chama constantemente sin() e cos() com um parâmetro que aumenta em etapas constantes (neste caso, $$$\ frac {2 \ pi} {\ text {num}}$) Este artigo é sobre como otimizar esse código para velocidade de forma que os mesmos cálculos possam ser executados sem o uso das funções sin() ou cos() (no loop interno) e até mesmo da função sincos() mais rápida.

Mas se você olhar para a primeira fórmula do artigo, veremos que, se $$$f = \ alpha$e $$$t = \ beta$podemos reescrever isso como

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

ou em outras palavras

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

Como f é uma constante, sin(f) e cos(f) também. Chame-os de b respectivamente:

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

Essa expressão pode ser usada diretamente em nosso código e, em seguida, obtemos um loop no qual funções trigonométricas caras (e de fato não) são chamadas!

 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; } 

Interpretação

Até agora, brincamos cegamente com matemática, sem entender o que realmente está acontecendo. Vamos reescrever o loop interno assim:

$$

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

Alguns de vocês devem ter notado que essa é uma fórmula para rotacionar um objeto no espaço bidimensional. Se você ainda não entende isso, talvez o formulário da matriz o ajude.

$$

$$\ 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)$$

De fato, sin(t) e cos(t) podem ser agrupados em um vetor de comprimento 1 e desenhados como um ponto na tela. Chame esse vetor x . Então $$$x = \ {\ sin \ beta, \ cos \ beta \}$ . Portanto, a forma vetorial de expressão é $$$x_ {n + 1} = Rx_n$onde $$$R = \ left (\ begin {matrix} a & b \\ - b & a \ end {matrix} \ right)$ . Vemos que nosso loop executa uma pequena rotação bidimensional a cada iteração, de modo que x gira em círculo durante a execução do ciclo. Cada iteração x gira em $$$\ frac {2 \ pi} {\ text {num}}$radiano.
Então basicamente

$$

$$\ 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$$

esta é a fórmula do movimento pontual $$$x = \ {\ cos \ alpha, \ sin \ alpha \}$ circunferência em incrementos de $$$\ beta$radiano. Para isso, construiremos um dos dois eixos de rotação, por exemplo, $$$u = \ {\ cos \ beta, \ sin \ beta \}$ . O primeiro componente da rotação é a projeção. $$$x$em $$$u$. Desde $$$x$e $$$u$normalizada (tem um comprimento de 1), a projeção é o produto escalar. Portanto, $$$s = x \ cdot u = \ sen \ alfa \ cdot \ cos \ beta + \ cos \ alfa \ cdot \ sin \ beta$e, é claro, o segundo componente é a anti-projeção, que pode ser encontrada projetando-se no eixo perpendicular, $$$v$. Podemos criar esse vetor expandindo as coordenadas $$$u$e mude o sinal para o oposto na primeira coordenada: $$$c = x \ cdot v = \ cos \ alpha \ cdot \ cos \ beta + \ sin \ alfa \ cdot \ sin \ beta$

Anotações

Normalmente, você deve poder executar essas pequenas rotações repetidas vezes. De fato $$$| R | = \ left | \ begin {matriz} a & b \\ - b & a \ end {matriz} \ right | = a ^ 2 + b ^ 2 = \ sen ^ 2 \ alfa + \ cos ^ 2 \ alpha = 1$ , o que significa que a matriz $$$R$não aumenta nem diminui o espaço no qual é aplicado, o que significa que $$$x$irá se mover em um círculo perfeito. No entanto, como os computadores não são precisos, $$$x$irá se mover em espiral e eventualmente coincidir com o centro do círculo de rotação. Eu não tive esses problemas, mas acho que eles podem ocorrer com um num muito grande, ou seja, pequenos ângulos de viragem.

Exemplo

No Kindercrasher , a demo de 4096 bytes de 2008 (captura de tela no KDPV), um grupo de esferas pulsa na música. Para isso, contei a transformação de Fourier do som. Existem algoritmos que fazem isso em tempo real, por exemplo, FFT . No entanto, meu código deve caber em alguns kilobytes e eu decidi seguir o outro caminho. Em vez de implementar a FFT, escrevi a DFT por sua definição simples. Você pode conferir na wikipedia.

$$

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

Minha função também utiliza um buffer de áudio estéreo de 16 bits, x , e retorna as primeiras 128 frequências do espectro sonoro do som y . Veja como o loop interno é organizado, o que executa 4096 vezes: nem uma única chamada para as funções sin() ou cos() , embora em outras implementações essas chamadas sejam.

 #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); } }