Inicia verificando se o ANGBUSCA é maior ou menor que ANGENTRADA. Inicialmente ANGBUSCA é zero, portanto será menor.
Quando MENOR, ele SOMARÁ o angulo apontado na tabela na variável ANGE, e também somará o XS em Y, e subtrairá YS de X.
Quando ANGBUSCA for MAIOR que ANGE, a coisa é ao contrário, ele SUBTRAIRÁ o angulo apontado na tabela na variável ANGE, agora SUBTRAIRÁ XS de Y, e SOMARÁ YS em X.
Se você entendeu isso, entendeu CORDIC.
Ele fará de 14 a 21 loops, de acordo com o que você quer de precisão e elementos na tabela.
A cada LOOP ele incrementa o ponteiro na tabela, apontando para um ângulo menor e menor, e usa esse ângulo menor para ser somado ou subtraido da variável ANGBUSCA, a fim de fazer um tipo de perseguição, buscando fazer essa variável ANGBUSCA ter o conteúdo mais próximo possível da ANGENTRADA. Como os valores na tabela são o resultado de uma divisão binária da Tangente, é óbvio que com o mesmo número de divisões na busca, a variável ANGEBUSCA terá um valor muuuito próximo da ANGENTRADA, com diferença muito pequena.
Também, a cada loop, conforme ANGEBUSCA seja menor ou maior, terá o valorzinho da tabela sendo somado ou subtraído à essa variável, e o valor de XS e YS será somado ou subtraído de Y e X.
Mas veja que X iniciou com o CORDIC Gain 0.60725... e isso é a multiplicação de todos os cosenos dos ângulos da tabela. Como o valor atual de X e Y em cada passada é shiftado para a direita, tantos bits quanto for o numero da passada, esse valor de X e Y representam o valor do coseno do angulo cujo valor da tangente foi também dividida binariamente. O que ele faz ao subtrair ou adicionar o XS e YX em Y e X (veja que é cruzado), é exatamente mantendo a relação entre seno e coseno também buscando ficar o mais próximo possível de ANGBUSCA. Quando seno sobe, coseno desce e vice-versa.
A mágica é muito boa. Confuso? olhe o código ao final e entenderá.
Ao completar as "n" passadas, o valor em X e Y corresponderão ao coseno e seno do angulo na variável ANGENTRADA, com resolução e precisão. Rodando 14 passadas já é suficiente para ter uma boa precisão, mas quanto mais passadas melhor. A quantidade de passadas também depende da escala do CORDIC Gain, usando 0.6072 x 65536 (16 bits) esse valor vai chegar a zero já após 14 ou 15 passadas, não adianta continuar fazendo passadas, com esse valor já zerado (XS e YS), nada será adicionado ou subtraído das variáveis X e Y, então para ter mais passadas, o CORDIC Gain terá que usar mais de 16 bits, para pelo menos 8 dígitos decimais de precisão será necessário usar 3 bytes (24 bits) no CORDIC Gain (e inicial em X).
Notar que, a quantidade de bits (multiplicador do CORDIC Gain) nada tem a ver com o multiplicador do angulo de entrada, nem da tabela de ângulos, esses dois tem que usar o mesmo multiplicador, mas 14 bits é suficiente (multiplicar o angulo de entrada e a tabela por 16384, por exemplo).
Mas é claro que, num microcontrolador você não trabalha com números tão pequenos, querendo trabalhar com números inteiros. Então tudo, desde a tabela, quanto ao CORDIC Gain, você já multiplica por 65536 para ter dois bytes inteiros a brincar. Eu já usei 3 bytes e multipliquei logo por 2^24.
Parte lógica do código, todas as variáveis são integers e unsigned, mas ocorre um probleminha de valor negativo se a ENGENTRADA for menor de 4.8... ou maior de 85.2... (algo assim), mas é facilmente resolvível usando 3 bytes nas variáveis X e Y e quando shiftar XS e YS abaixo, o primeiro shift ser aritmético, ou seja, copiar o bit de alta ordem para os bits sendo inicialmente empurrados para a direita, assim, se o valor de X ou Y ficarem negativos, o valor dividido por SHIFTS (shift right) também será negativo, e a soma ou subtração se encarrega de resolver esse negativo mesmo usando números inteiros unsigned. Assembly é lindo. O mundo todo usa variáveis signed para o CORDIC, exatamente por esse problema dos ângulos pequenos, eu resolvi isso somando 5 graus no ANGENTRADA e também inserindo 5 * 65536 em X no início, e usando o artifício do aritmetic shift right em XS e YS mantendo o bit de mais alta ordem significando o sinal shiftado, pronto, resolvido. Se X = FF0224, o bit de alta ordem me diz que ele é negativo, pois nenhum valor positivo chegará a ser tão alto em 3 bytes (por mais que seja 90°x65536 será 5A0000, bit de alta ordem será sempre zero para positivos). Ao usar ASR para shiftar para a direita, ficará FF8112, mantendo a informação negativa, e ao somar, ou subtrair (de Y) resolverá sozinho. Se não usasse ASR, após o shift seria 7F8112, o que seria considerado um valor positivo na soma ou subtração e causaria problemas.
Abaixo o bicho gigante, cabeludo, cheio de dentes e unhas, o famoso CORDIC para calcular seno e coseno.
Ah sim, usando o Mult = 65536, o valor resultando de seno e coseno em Y e X também estará entre 0 e 65536, ou seja, é só considerar o ponto decimal em frente ao número obtido. Se ANGENTRADA for 30, a variável Y terminará com o valor 32768, que significa 0.5 (metade) de 65536.
Depois dê uma olhadinha nos outros websites da Internet explicando CORDIC e veja como lá o bicho é pintado de cabeludo, como por exemplo esses:
Código simples em Basic para entender:
Mult = 65535
Tab(1) = 45 * Mult
Tab(2) = 26.56505 * Mult
Tab(3) = 14.03624 * Mult
Tab(4) = 7.125016 * Mult
Tab(5) = 3.576334 * Mult
Tab(6) = 1.789911 * Mult
Tab(7) = 0.895174 * Mult
Tab(8) = 0.4476142 * Mult
Tab(9) = 0.2238105 * Mult
Tab(10) = 0.1119057 * Mult
Tab(11) = 0.0559529 * Mult
Tab(12) = 0.0279765 * Mult
Tab(13) = 0.0139822 * Mult
Tab(14) = 0.0069941 * Mult
Tab(15) = 0.0034971 * Mult
Tab(16) = 0.0017485 * Mult
Tab(17) = 0.0008743 * Mult
Tab(18) = 0.0004371 * Mult
Tab(19) = 0.0002186 * Mult
Tab(20) = 0.0001093 * Mult
X = 0.60725293 * Mult
Y = 0
ANGBUSCA = 0
ANGENTRADA = 30 ; entrada do usuário
XS = 0
YS = 0
SHIFTS = 0
LOOP:
FOR L = 1 to 20
YS = Y / SHIFTS
XS = X / SHIFTS
IF ANGENTRADA > ANGBUSCA THEN
X = X - YE
Y = Y + XE
ANGBUSCA = ANGBUSCA + TAB(L)
ELSE
X = X + YE
Y = Y - XE
ANGBUSCA = ANGBUSCA - TAB(L)
END IF
SHIFTS = SHIFTS * 2
NEXT L
Vamos fazer um exercício prático, e multiplicando tudo por 65536 para ter inteiros.
Angulo de entrada = 30° * 65536 = 1966080
CORDIC Gain = 0.607252... x 65536 = 39797
Angulos somados ou subtraídos de AngBusca abaixo, são os valores da tabela x 65536.
Novamente, na lista abaixo, a cada passo XS = X/2^passo e YS = Y/2^passo
Vermelho AngBusca = AngBusca + Tabela, Y = Y + Xs, X = X - YS
Azul AngBusca = AngBusca - Tabela, Y = Y - XS, X = X + YS
AngEntrada AngBusca X Y XS YS
1966080 0 39797 0 39797 0 Passo 0
1966080 2949120 39797 39797 19899 19899 Passo 1
1966080 1208153 59696 19898 14924 4975 Passo 2
1966080 2128032 54721 34822 6840 4353 Passo 3
1966080 1661087 59074 27982 3692 1749 Passo 4
1966080 1895466 57325 31674 1792 990 Passo 5
1966080 2012770 56335 33466 880 523 Passo 6
1966080 1954104 56858 32586 444 255 Passo 7
1966080 1983439 56603 33030 221 129 Passo 8
1966080 1968771 56732 32809 111 64 Passo 9
1966080 1951437 56796 32698 56 32 Passo 10
1966080 1965104 56764 32754 28 16 Passo 11
1966080 1966937 56748 32782 14 8 Passo 12
1966080 1966020 56756 32768 7 4 Passo 13
1966080 1966478 56752 32775 4 2 Passo 14
1966080 1966249 56754 32771 2 1 Passo 15
1966080 1966134 56755 32769 1 1 Passo 16
1966080 1966077 56756 32768 0 0 Passo 17
Valor Final Y = Sin(30) = 32768 / 65536 = 0.5
Valor Final X = Cos(30) = 56756 / 65536 = 0.866027...
Coincidentemente os valores corretos já foram obtidos no passo 13, mas é impossível dizer que eram corretos antes de usar toda a tabela e também o ANGBUSCA estava com 40 de diferença do ANGENTRADA.
Esse mecanismo do CORDIC pode ser usado (com outros valores de tabela e Gain) para calcular logaritmos, e outros. Pense como você faria a tabela para calcular números elevados a potência, ou raiz quadrada.
Fato também interessante, é que o algoritmo CORDIC acima, tanto funciona em Graus como em Radianos, basta que a tabela seja gerada na escolha, e o angulo de entrada também.
Mais tarde posto aqui o código assembly do AVR para calcular sin(x) cos(x).
OTIMIZANDO O ALGORITMO CORDIC ORIGINAL
Preparando e já otimizando o CORDIC original, para implementar em AVR Assembly, existem alguns detalhes que já podemos ir mudando, para acelerar o processo. Inicialmente vejo 3 possíveis melhorias.
OTIMIZAÇÃO #1
Por exemplo, observe que o ANGBUSCA inicia com zero, e obviamente ele será SEMPRE menor que o ANGENTRADA no primeiro passo. Então, por que não eliminar o "Passo 0" e "Passo 1", se já sabemos o que irá acontecer? Essa alteração economizaria em torno de 5~6% do tempo total.
Veja a linha do "Passo 1", ANGBUSCA vai receber 45°, o primeiro ítem da Tabela, aqui no caso multiplicado por 65536 = 2949120. Também, veja que quando ANGBUSCA é menor que ANGENTRADA e soma valor da tabela no ANGBUSCA, também soma o valor de X em Y e já prepara XS e YS com a metade de X e de Y, para somar ou subtrair em Y e X, conforme se ANGBUSCA ficar maior ou menor que ANGENTRADA no "Passo 2".
Então, o normal é X receber o valor de CORDIG GAIN, aqui no caso 39797 e Y receberia zero, mas como estamos eliminando o "Passo 0", já colocamos CORDIG GAIN também em Y.
Então, para eliminar o "Passo 0", faremos X = Y = 39797, ANGBUSCA = 2949120, XS = YS = 19899, e já partimos para o "Passo 2".
OTIMIZAÇÃO #2
Outra melhoria a ser feita, que funciona bem em Assembly, é evitar comparar ANGENTRADA com ANGBUSCA a cada passo, isso requer 4 instruções, 3 para comparar os registradores que compõe ANGBUSCA com ANGENTRADA e a instrução de salto condicional. Podemos então no início carregar as variáveis ANGENTRADA com zero e ANGBUSCA com o valor do ângulo de entrada, porém negativo. Ou seja, se o angulo a calcular o sin() for 30°, vezes 65536 = 1966080 (0x1E0000) o negativo é 0xE20000, a carregar no ANGBUSCA. Motivo dessa mudança? Ao somar os valores da Tabela a cada passada, se o resultado da soma for valor positivo, significa que ANGBUSCA acabou maior que ANGENTRADA (que é zero), se ficar negativo significa que ficou MENOR que ANGENTRADA. Então, não precisa comparar os 3 bytes e decidir, basta olhar o bit de mais alta ordem do byte de mais alta ordem do ANGBUSCA (no caso acima, 0xE2 = 0b11110010, se estiver levantado, é porque o valor é negativo, então ANGBUSCA é menor que ANGENTRADA, e vice-versa. Também podemos observar o C (carry) bit da soma ou subtração do byte de mais alta ordem (igual ao bit 7 desse byte) A decisão seria só por verificar um bit. No AVR existe a instrução Skip if Bit in Register is Cleared/Set (SBRC, SBRS), seria só uma instrução ao invés de quatro, vezes a quantidade de passos. Se usarmos 21 passos, seriam 63 instruções de economia, se rodando a 8MHz economizaríamos 7.8 microsegundos, isso é quase 6% to tempo total.
Ou seja, nesse jogo de somar e subtrair valores da Tabela no registrador ANGBUSCA, o programa estaria tentando manter ANGBUSCA mais próximo possível de zero.
Você pode se perguntar, qual é a grande vantagem de fazer tudo isso só para economizar 15~20 microsegundos? Veja, quando escrevendo código em assembly, você busca excelência, rapidez, menor tamanho, facilidade de compreensão da rotina, extremamente funcional, preciso e livre de erros. Se não estiver preocupado com isso, escreva em linguagem de alto nível, como C, e use bibliotecas prontas, e aceite o resultado que for.
OTIMIZAÇÃO #3
O algoritmo tradicional CORDIC compara ANGENTRADA com ANGBUSCA e decide se soma ou subtrai o próximo valor da tabela em ANGBUSCA, e se soma e subtrai valores XS e YS de Y e X, veja abaixo:
SHIFTS = 0
FOR L = 1 to 20
YS = Y / SHIFTS
XS = X / SHIFTS
IF ANGENTRADA > ANGBUSCA THEN
X = X - YE
Y = Y + XE
ANGBUSCA = ANGBUSCA + TAB(L)
ELSE
X = X + YE
Y = Y - XE
ANGBUSCA = ANGBUSCA - TAB(L)
END IF
SHIFTS = SHIFTS * 2
NEXT L
Também observe acima que a cada passo, XS e YS é o resultado de X e Y deslocado para a direita 2^(L-1) vezes, a variável SHIFTS acima representa isso. Veja na lista de resultados, passos 8 a 12 abaixo, observe passo 10, X = 56796 (0xDDDC) será deslocado 10 bits para a direita para gerar XS = 56 (0x38). Também veja que no espaço entre XS = 55.5 a 56.49, cujo inteiro é 56, ao elevar esses décimos 2^10, ou exatamente a variação 1 x 2^10, teremos uma possível variação em X de 1024. Isso significa que não importa se o valor 64 de YS (YE no loop acima) no passo 9 seja somado ou subtraído de X, que não vai mudar o valor em XS. Ou seja, é garantido que, do passo 10 até o final o valor exato de YS somado ou subtraído de X não irá interferir mais em XS, então a partir desse passo basta deslocar o valor anterior de XS e YS um bit à direita a cada passo.
O resumo disso é que a partir do passo 10 (metade dos passos) já não é mais necessário deslocar X e Y (passos-1) à direita para obter YS e XS, será só perda de tempo e ciclos de clock. A maneira fácil de verificar qual passo é o limite, é verificar se 2^(passo-1) é maior ou menor que XS na tabela de resultados. Como no exemplo foram usados 16 ou 17 passos, veja que 2^8 = 256 e é maior que 221 (XS no passo 8), mas 2^7 (128) é menor que XS 444 no passo 7, portanto já a partir do passo 8 já faz sentido aplicar só um deslocamento à direita dos XS e YS anteriores para os novos XS e YS e economizar ciclos de clock. Perceba que o valor 8 é metade do numero de passos. Se usando 22 passos, já a partir do passo 11 é possível usar essa otimização.
Como a quantidade de deslocamentos dos bits de XS e YS variam de acordo com o número do passo, e como essa otimização é usada da metade mais alta de passos, estaríamos economizando a maior parte dos deslocamentos. Como XS e YS no exemplo usam 16 bits, são dois deslocamentos e o loop do número de passos, então, do passo 8 até o 16 seriam 8+9+10+11+12+13+14+15+16 = 108 deslocamentos em 2 bytes + um contador + jump relativo a cada passo, a transferência de X para XS e Y para YS antes dos deslocamentos, a economia seria em torno de 6 clock cycles por deslocamento, em torno de 650 clock cycles, em 8MHz significa ~81 microsegundos. Porém agora teríamos deslocamento simples de 2 bytes x 8 = 16 ciclos de clock adicionais, economia final em torno de ~78 microsegundos. Essa economia significa em torno de 48% do tempo total, somado com os 5 e 6% das otimizações anteriores, estaríamos chegando a uma redução de ~58% da rotina original do CORDIC.
Em assembly, nos passos acima de 8 já podemos simplesmente mover Byte2 para Byte3 e Byte1 para Byte2, ao invés de fazer 8 deslocamentos entre os 3 bytes, isso já reduz um bom percentual de ciclos de clock no processo, mas a otimização #3 reduz bem mais.
X Y XS YS
1966080 1954104 56858 32586 444 255 Passo 7
1966080 1983439 56603 33030 221 129 Passo 8
1966080 1968771 56732 32809 111 64 Passo 9
1966080 1951437 56796 32698 56 32 Passo 10
1966080 1965104 56764 32754 28 16 Passo 11
1966080 1966937 56748 32782 14 8 Passo 12
Resumo das otimizações 1, 2 e 3:
Eliminar passo 0 e 1, carregando os registradores com valores já conhecidos a acontecer.
Mudar a forma de comparar se o ANGBUSCA é maior ou menor que ANGENTRADA, fazendo ANGBUSCA = ANGENTRADA (negativo) e zerando ANGENTRADA, e só observando o bit de alta ordem de ANGBUSCA para decidir se usa as instruções em azul ou vermelho da rotina.
IF ANGENTRADA > ANGBUSCA THEN
X = X - YE
Y = Y + XE
ANGBUSCA = ANGBUSCA + TAB(L)
ELSE
X = X + YE
Y = Y - XE
ANGBUSCA = ANGBUSCA - TAB(L)
END IF
E por fim, eliminar os deslocamentos de X e Y para XS e YS 2^(passo-1) a partir da metade dos passos até o fim, fazendo só um simples shift-right do valor anterior de YS e XS.
.jpg/640px-Capacitors_(7189597135).jpg)
