CORDIC - ESSE TIGRE NÃO MORDE
CORDIC é um procedimento usado para determinar e calcular valores como seno(x), coseno(x), logarítmo(x), tangente(x) e outros, sem usar multiplicação ou divisão, nem formulas malucas ou tabelas gigantes.
CORDIC é usado em muitas calculadoras e até programas de computador, onde velocidade é importante, mas sem requerer muito hardware para tais cálculos.
A criação do CORDIC teve por base exatamente isso, velocidade e facilidade, e acabou sendo tão simples, que até mesmo soluções mecânicas foram inventadas para tal.
Pois é, depois de anos com medo desse bicho, fui atrás e cutuquei com vara curta.
Pensei comigo mesmo que já era hora de ao menos entender como funciona nos detalhes.
Ao ler os textos na Internet, me apavorei ainda mais, o povo usa formulas complexas, e desenhos que confundem mais do que explicam. Mas não desisti, pois sempre pensei, "como pode algo simples ser tão complexo e com descrição tão assustadora?".
Devagar e ao meu tempo fui entendendo e montando um racional mais simples de explicar como o CORDIC funciona, de forma a poder escrever um programa para o microcontrolador AVR em assembly. Para quem conhece assembly, é claro que você precisa entender muito bem e em detalhes o que está fazendo, portanto, quando eu decidi escrever CORDIC para o AVR em assembly, é porque eu já tinha entendido tal procedimento muito bem. Descobri que CORDIC não morde, na verdade é até manso demais.
No passado escrevi uma rotina, também em AVR assembly, para calcular o valor de seno(x) e coseno(x), considerando o ângulo com uma casa decimal, de zero a 90 graus (0 ~ 89.9°), usando o sistema de aproximação de Taylor que usa uma pancada de fator de potência e divisões por fatoriais. Foi complicado e demorou. Agora gostaria de comparar a performance daquela rotina contra o CORDIC, e por isso meti as caras no programa.
Primeiro, CORDIC faz uma mágica com uma tabela pequena, de no máximo 22 ou 23 entradas. Essa tabela contém os valores de ângulos entre zero e 45°. Inicia em 45° e vai decrementando até chegar a valor muito pequeno, a última entrada é 0.00005464°. Quase tudo em linguagem de programação usa divisão binária, por dois, mas não nessa tabela. Abaixo os primeiros 10 elementos da tabela;
45.000000
26.565051
14.036243
7.125016
3.576334
1.789911
0.895174
0.447614
0.223811
0.111906
...
Esses valores de ângulos são calculados antecipadamente no Excel ou calculadora e digitados no programa, eles representam o ângulo cuja tangente é 1, e dividida por dois consecutivamente. Ou seja, a atan(1/2^n) onde n vai de zero a o valor máximo de elementos na tabela. Os valores foram calculados com base em:
atan(1/1) ==> atan(1)
atan(1/2) atan(0.5)
atan(1/4) atan(0.25)
atan(1/8) atan(0.125)
atan(1/16) atan(0.0625)
...
Note que a divisão por dois não é do ângulo, mas da tangente, e a tabela mostra os angulos que representam essa tangente sendo divida por dois. Isso é fundamental para entender CORIC. Então a tabela contém os ângulos que representam a divisão binária do valor da tangente entre 45° e praticamente zero.
Mas o que a tangente tem a ver com isso?
Lembrar algumas propriedades de trigonometria:
Sin(x)^2 + Cos(x)^2 = 1
Por exemplo, Sin(30) = 0.5, Cos(30) = 0.866, 0.5^2 = 0.25, 0.866^2 = 0.75, somando ambos = 1.
O mais importante aqui, a Tangente = Seno(x) / Coseno(x).
Então, sabendo a tangente e o ângulo, já é meio caminho andado para saber o seno(x) e coseno(x).
Mas como saber a tangente se só se tem o ângulo?
Aha, ai é que entra a tabela.
Mas a tabela só tem poucos valores de ângulo que corresponde à tangente.
Ai é que entra a mágica o número K, que é conhecido como "ganho do CORDIC".
Multiplicando todos os cosenos dos ângulos da tabela, ou seja, cos(45) * cos(26,56505) * cos(14.03624)... ter/a
um valor de 0.6072529350088810... e quanto mais preciso esse numero, mais precisos serão os valores obtidos de seno(x) e coseno(x) pelo algoritmo CORDIC. Esse valor ai acima será sempre o mesmo, nem precisa mais calcular, pode usa-lo sem medo. CORDIC Gain. Na internet você fica louco tentando entender o que significa esse CORDIC Gain, o povo complica demais a explicação. Só como fato curioso, o ângulo cujo coseno é 0.6072... é 52.60886...°, e a princípio não tem nada a ver com coisa alguma, mas é isso mesmo.
Já estamos chegando na mágica. Veja que existe uma relação entre os valores da tabela, que são os ângulos correspondentes a tangente de 45° que é "1.00", e tangentes sendo divididas binariamente por dois, e agora os cosenos dos ângulos da tabela multiplicados um pelo outro, obtendo 0.6072... Estamos preparando uma divisão binária entre ângulos, tangentes e cosenos. A mágica.
O algoritmo CORDIC inicia com duas variáveis, uma de seno e outra de coseno, que chamaremos de Y e X. A variável de seno (Y) é iniciada com zero. A do coseno (X) é iniciada com o valor 0.6072... Também teremos duas outras variáveis, a do ângulo de entrada (o que o usuário entrou para buscar o seno e coseno), e a do ângulo-de-busca. Vamos dar nomes, a do ângulo de entrada, chamaremos de ANGENTRADA, e a do ângulo de busca, de ANGBUSCA. A ANGENTRADA será iniciada com o ângulo do usuário, a ANGBUSCA com zero. Também um ponteiro apontando para o início da tabela, onde terá 45° como primeiro elemento. Outras duas variáveis de nome YS e XS, que recebem os valores de Y e X (atuais) shiftados à direita (divididos por dois), a cada rodada da rotina
Pronto, o CORDIC está pronto para rodar e é muito simples.
O que ele faz é um joguinho de ping-pong e uma busca binária.
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:
ou até mesmo na Wikipedia (que deveria ser simples): https://en.wikipedia.org/wiki/CORDIC
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.
Nenhum comentário:
Postar um comentário