{twitter}: @FpgaWarrior

CORDIC es un acrónimo de COordinate Rotation DIgital Computer que en castellano significa Computadora Digital para Rotación de Coordenadas. El algoritmo original fue propuesto por Jack Volder en el año 1959, con el propósito de calcular funciones trigonométricas mediante la rotación de vectores. La rotación de vectores puede utilizarse a su vez para la conversión de sistemas de coordenadas (cartesiano a polar y viceversa), magnitud de vectores y como parte de funciones matemáticas mas complejas como pueden ser la Transformada Rápida de Fourier (Fast Fourier Transform o FFT) y la Transformada Discreta del Coseno (Discrete Cosine Transform o DCT). El algoritmo CORDIC, proporciona un método iterativo para rotar vectores en ciertos ángulos determinados y se basa exclusivamente en sumas y desplazamientos. El algoritmo original tal como fue propuesto inicialmente, describe la rotación de un vector bidimensional en el plano cartesiano. Su funcionamiento se deduce de la fórmula general para rotación de vectores como se muestra a continuación:

La siguiente fórmula se deduce a partir de las identidades trigonométricas conocidas. Si r = 1 y θ = ω −φ entonces ω = θ + φ. Utilizando las identidades para el seno y el coseno de la suma de ángulos se tiene:

Si se sustituye ω = θ + φ se puede obtener:

Reagrupando los términos anteriores y considerando que cosθ ≠ 0:

En el algoritmo CORDIC, las rotaciones de vectores son reemplazadas por pseudo-rotaciones como se muestra en la figura siguiente:

El algoritmo CORDIC se opera normalmente en dos modos. El primero, se denomina rotación (rotation) y consiste en rotar el vector de entrada un ángulo específico que se introduce como parámetro. El segundo modo, denominado vectorización (vectoring), rota el vector de entrada hacia el eje X, acumulando el ángulo necesario para efectuar dicha rotación. En el caso de una rotación, el acumulador angular se inicializa con el ángulo a rotar. La decisión sobre el sentido de rotación en cada paso de iteración, se efectúa para minimizar la magnitud del ángulo acumulado. Por ello, el signo que determina el sentido de rotación, se obtiene del valor de dicho ángulo en cada paso. Para el modo rotación, las ecuaciones son:

Para una vectorización, el ángulo ingresado se rota para alinearlo con el eje X. Para obtener este resultado, en lugar de minimizar la magnitud del acumulador angular, se minimiza la magnitud del componente y, ya que si y = 0 entonces el vector se encuentra sobre el eje X. Asimismo se utiliza el signo del componente y para determinar la dirección de rotación. Si el acumulador angular se inicializa con cero, al final del proceso contendrá el ángulo de rotación adecuado. Por lo tanto se pueden deducir las siguientes ecuaciones:

Durante algún tiempo hubo una tendencia que buscaba implementar el cálculo de las funciones trigonométricas sin() y cos() empleando series de Taylor. El problema con este enfoque es que requiere de la realización de muchas multiplicaciones. CORDIC tiene un rendimiento superior al enfoque anterior y es muy utilizado en implementaciones de hardware cuando el sistema en cuestión no dispone de un módulo para realizar multiplicaciones (ya sea un micro-controlador o una FPGA) o cuando el área efectiva de chip en una FPGA no dispone de espacio suficiente para implementar un módulo de multiplicación. Además existen varias arquitecturas a seguir a la hora de implementar el algoritmo CORDIC.

Arquitectura Bit-Paralela Iterativa:

Cuando el algoritmo CORDIC se implementa en hardware, una de las arquitectura posibles es la arquitectura bit-paralela iterativa. La denominación de paralela se debe a la forma en que se opera con las componentes X, Y y Z. En esta arquitectura cada etapa del algoritmo consiste de un registro para almacenar la salida, una unidad de desplazamiento y un sumador algebraico. Al comenzar el cálculo, los valores iniciales para 0 x , 0 y y 0 z ingresan en forma paralela a los registros a través del multiplexor. El bit más significativo del componente Z ó Y en cada paso de iteración determina la operación a efectuar por el sumador algebraico en modo rotación o vectorización respectivamente. Las señales correspondientes a los componentes X e Y son desplazadas y luego restadas o sumadas a las señales sin desplazar, correspondientes al componente opuesto. El componente Z combina los valores almacenados en el registro con valores que obtiene de una tabla de búsqueda (Lookup Table, LUT) de arcotangentes precalculadas con una cantidad de entradas proporcional a la cantidad de iteraciones. Esta tabla de búsqueda puede ser implementada como una memoria ROM. La dirección de memoria cambia acorde al número de iteración. Para n iteraciones el valor a calcular puede ser obtenido en la salida. Se requiere de un controlador que puede ser implementado como una máquina de estados para controlar los multiplexores, la cantidad de desplazamientos y el direccionamiento de las constantes precalculadas. Esta representación tiene como ventaja el uso eficiente de hardware, debido a que los recursos (registros, multiplexores, unidades de desplazamiento y sumadores) son reutilizados en cada iteración. En la siguiente figura se muestra el esquema correspondiente a la arquitectura descripta. La señal Modo especifica el modo de operación (Rotación o Vectorización). Se considera además en el esquema que las componentes tienen un ancho de m bits. Sin embargo se puede ampliar el ancho en bits de los buses internos para obtener una mejor exactitud.

Arquitectura Bit-Paralela Desplegada:

En lugar de almacenar el resultado de cada paso de iteración en registros y volver a utilizar los mismos recursos, el diseño puede desplegarse como muestra siguiente figura. El diseño se separa en etapas correspondientes a cada iteración. Cada etapa está compuesta por los mismos componentes, dos unidades de desplazamiento y dos sumadores algebraicos. Por consiguiente la salida de una etapa corresponde a la entrada de la siguiente etapa. Los valores iniciales para X0 , Y0 y Z0 se entran en paralelo a la primera etapa. En el esquema se supone un ancho de palabra de m bits. La señal Modo indica el modo de operación al igual que en el caso iterativo.

Este diseño introduce dos ventajas importantes. Las unidades de desplazamiento así como las constantes correspondientes a cada iteración pueden ser cableadas. La segunda ventaja radica en que el circuito es puramente combinatorio, y por lo tanto no se necesita de una unidad de control, simplificando enormemente el diseño. La desventaja que presenta esta arquitectura es la enorme cantidad de espacio que requiere su implementación.

Arquitectura Bit-Serie Iterativa:

En un diseño en serie, se procesa un bit por vez, con lo cual las interconexiones se reducen al ancho de un bit. El desempeño δ de esta arquitectura se puede calcular aplicando la siguiente expresión:

La arquitectura resultante se muestra en la siguiente figura. Cada etapa del algoritmo consiste en un multiplexor, un registro de desplazamiento y un sumador algebraico bit-serie. El sumador algebraico se implementa como un sumador completo en el que una resta se lleva a cabo sumando el complemento a dos del valor a restar. La operación a efectuar por el sumador algebraico se obtiene del bit de signo del componente z para el modo rotación. Para el modo vectorización, del bit de signo se obtiene del componente y.

La operación de desplazamiento (multiplicación por 2^-i ) se lleva a cabo leyendo el bit i – 1 considerado a partir del extremo derecho del registro de desplazamiento. Se puede usar un multiplexor para cambiar la posición de acuerdo a la iteración actual. Los valores iniciales para 0 x , 0 y y 0 z ingresan al registro de desplazamiento por el extremo izquierdo en el esquema. Cuando el primer bit que fue ingresado al registro es procesado por el sumador algebraico, el multiplexor permite nuevamente el ingreso de los bits sumados al registro de desplazamiento. Por último, cuando se han completado todas las iteraciones los multiplexores permiten el ingreso de un nuevo valor al registro de desplazamiento y el valor calculado se obtiene en la salida. Alternativamente la carga y lectura del registro de desplazamiento puede efectuarse en paralelo para simplificar el diseño. La desventaja que presenta esta arquitectura con respecto a las que procesan los vectores de entrada en forma paralela es que se introduce un retardo proporcional al ancho de palabra en cada etapa, ocasionado por los desplazamientos. Por otra parte se requiere de hardware de control mas complejo. La ventaja que presenta esta arquitectura es que los buses de interconexión entre los componentes son del ancho de un bit y el registro de desplazamiento está integrado con el registro intermedio, minimizando espacio al momento de su implementación.

Para finalizar esta entrada voy a poner a continuación algunos ejemplos de código que pueden servir para adentrarse en formas o métodos de implementación de este algoritmo. El siguiente bloque de código es una implementación de Matlab del algoritmo CORDIC que no depende de ninguna función trascendental excepto en el cálculo previo de tablas (LookUp Table). Si el número de iteraciones n está predeterminado, entonces la segunda tabla puede reemplazarse por una única constante.

function v = cordic(beta,n)
% This function computes v = [cos(beta), sin(beta)] (beta in radians)

if beta < -pi/2 || beta > pi/2
    if beta < 0
        v = cordic(beta + pi, n);
    else
        v = cordic(beta - pi, n);
    end
    v = -v; % flip the sign for second or third quadrant
    return
end

angles =  [  ...
    0.78539816339745   0.46364760900081   0.24497866312686   0.12435499454676 ...
    0.06241880999596   0.03123983343027   0.01562372862048   0.00781234106010 ...
    0.00390623013197   0.00195312251648   0.00097656218956   0.00048828121119 ...
    0.00024414062015   0.00012207031189   0.00006103515617   0.00003051757812 ...
    0.00001525878906   0.00000762939453   0.00000381469727   0.00000190734863 ...
    0.00000095367432   0.00000047683716   0.00000023841858   0.00000011920929 ...
    0.00000005960464   0.00000002980232   0.00000001490116   0.00000000745058 ];

Kvalues = [ ...
    0.70710678118655   0.63245553203368   0.61357199107790   0.60883391251775 ...
    0.60764825625617   0.60735177014130   0.60727764409353   0.60725911229889 ...
    0.60725447933256   0.60725332108988   0.60725303152913   0.60725295913894 ...
    0.60725294104140   0.60725293651701   0.60725293538591   0.60725293510314 ...
    0.60725293503245   0.60725293501477   0.60725293501035   0.60725293500925 ...
    0.60725293500897   0.60725293500890   0.60725293500889   0.60725293500888 ];
Kn = Kvalues(min(n, length(Kvalues)));

v = [1;0]; % start with 2-vector cosine and sine of zero
poweroftwo = 1;
angle = angles(1);

% Iterations
for j = 0:n-1;
    if beta < 0
        sigma = -1;
    else
        sigma = 1;
    end
    factor = sigma * poweroftwo;
    
    R = [1, -factor; factor, 1];
    v = R * v; % 2-by-2 matrix multiply
    beta = beta - sigma * angle; % update the remaining angle
    poweroftwo = poweroftwo / 2;
    % update the angle from table, or eventually by just dividing by two
    if j+2 > length(angles)
        angle = angle / 2;
    else
        angle = angles(j+2);
    end
end

% Adjust length of output vector to be [cos(beta), sin(beta)]:
v = v * Kn;
return

endfunction

A continuación se presenta una posible implementación del algoritmo CORDIC en lenguaje de programación JAVA que realiza 32 iteraciones para pre-calcular la tabla (LookUp Table). El algoritmo funciona con radianes, por lo que se debe volver a calcular el ángulo si el valor dado se proporciona en grados. A mayor número de iteraciones mejor precisión se obtiene.

public class Cordic {
 
    private double[] lookup;
    private final int iterations = 32;
    private final double K = 0.60725293510314;
    private final double halfPi = 1.57079632679;
 
    Cordic() {
        createLookupTable();
    }
 
    public double sin(double theta) {
        return sinCos(theta)[1];
    }
 
    public double cos(double theta) {
        return sinCos(theta)[0];
    }
 
    public double[] sinCos(double theta) {
        if (theta < 0 || theta > halfPi) {
            throw new IllegalArgumentException("Required value 0 < x < Pi/2");
        }
        double x = K;
        double y = 0;
        double z = theta;
        double v = 1.0;
        for (int i=0; i<iterations; i++) {
            double d = (z >= 0)? +1 : -1;
            double tx = x - d * y * v;
            double ty = y + d * x * v;
            double tz = z - d * lookup[i];
            x = tx; y= ty; z = tz;
            v *= 0.5;
        }
        double[] result = {x,y};
        return result;
    }
 
    private void createLookupTable() {
        lookup = new double[iterations];
        for (int i=0; i<iterations; i++) {
            lookup[i] = Math.atan(1 / Math.pow(2,i));
            System.out.format("Tan (%02d): %.14f / %018.3f %n", i,  lookup[i], Math.pow(2,i));
        }
    }
 
    public static void main(String[] args) {
        Cordic c = new Cordic();
        for (double i=0; i<90; i++) {
            double rad = i * (Math.PI / 180);
            System.out.format(
            "Sin: %04.1f (rad: %.14f) CORDIC: %.14f / java: %.14f %n", i, rad, c.sin(rad), Math.sin(rad) );
        }
    }
 
}

See you Keyboard CowBoy.