Algoritmus CORDIC s hodnotami uloženými ve formátu FX

Obsah

1. Algoritmus CORDIC s hodnotami uloženými ve formátu pevné řádové binární čárky
2. Datové typy, makra a pomocné funkce pro práci s formátem FX
3. Funkce určená pro naplnění tabulky atans[]
4. Výpočet tangenty pomocí algoritmu CORDIC – neoptimalizovaná verze
5. Výpočet tangenty pomocí algoritmu CORDIC – optimalizovaná verze
6. Výpis celého demonstračního příkladu
7. Obsah poslední části tohoto seriálu

1. Algoritmus CORDIC s hodnotami uloženými ve formátu pevné řádové binární čárky

V předchozích částech tohoto seriálu jsme si uvedli, jakým způsobem je možné algoritmus CORDIC (COordinate ROtation DIgital Computer) použít pro výpočet goniometrických funkcí s hodnotami uloženými ve formátu plovoucí řádové binární čárky (FP). Tentýž algoritmus je však po mírné modifikaci možné použít i při práci s formátem pevné řádové binární čárky (FX) a dá se říci, že teprve zde se plně ukazuje jeho implementační jednoduchost a současně i velká vyjadřovací síla – pomocí jednoho iteračního postupu je možné vypočítat poměrně velké množství navzájem značně odlišných funkcí, přičemž paměťové nároky algoritmu jsou minimální a rovnají se tabulce s cca deseti až dvaceti prvky s rozsahem 32 (popř. pouze 16) bitů. To je výhodné především při implementaci na mikrořadičích (oblíbená řada 8051, PICy, řada 68HC11 atd.) a také při práci s programovatelnými obvody typu FPGA.

2. Datové typy, makra a pomocné funkce pro práci s formátem FX

Před vlastním uvedením implementace algoritmu CORDIC s čísly uloženými ve formátu pevné řádové binární čárky je zapotřebí provést náležitou přípravu, protože formát FX (většinou) není ani mikroprocesory ani překladači podporován. Z tohoto důvodu si v této kapitole pro připomenutí uvedeme všechny základní aritmetické funkce určené pro práci s formátem FX. Jedná se o funkce fx_add(), fx_sub(), fx_mul() a fx_div() (tyto funkce byly podrobněji vysvětleny v předchozím díle). Kromě toho jsou uvedena i těla dalších pomocných funkcí, zejména fx_print() i maker určených pro převod stupňů na radiány a naopak. Pomocí konstant A a B je určen rozsah a přesnost numerických hodnot uložených ve formátu fx – ten je v naší implementaci shodný s třicetidvou­bitovým datovým typem signed int, jelikož potřebujeme pracovat jak s kladnými, tak i se zápornými hodnotami.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/* počet míst před a za binární řádovou tečkou */
#define A 16
#define B 16

/* Ludolfovo číslo */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* maximální počet iterací při běhu algoritmu */
#define MAX_ITER 16

/* "zesílení" při rotacích */
#define K 0.6073

/* převody mezi stupni a radiány */
#define rad2deg(rad) ((rad)*180.0/M_PI)
#define deg2rad(deg) ((deg)/180.0*M_PI)

/* datový typ, se kterým budeme pracovat */
typedef signed int fx;

/* hlavičky použitých funkcí */
void   fx_print(fx x);
fx     fp2fx(double x);
double fx2fp(fx x);

/* tabulka arkustangentu úhlů */
fx atans[MAX_ITER];

/* tabulka záporných celočíselných mocnin hodnoty 2 */
fx pows[MAX_ITER];

/*
 * Tisk numerické hodnoty uložené ve formátu pevné
 * řádové binární čárky (FX)
 */
void fx_print(fx x)
{
    int i;
    int val=x;                         /* pomocná proměnná pro převod do dvojkové soustavy */
    printf("bin: ");
    for (i=0; i<A+B; i++) {            /* převod na řetězec bitů (do dvojkové soustavy) */
        putchar(!!(val & (1<<(A+B-1)))+'0'); /* výpis hodnoty aktuálně nejvyššího bitu */
        if (i==B-1) putchar('.');      /* po řádové binární čárce vypsat značku */
        val=val<<1;                    /* posun na další (méně významný) bit */
    }

    printf("   hex: %08x   fp: %+11.5f\n", x, fx2fp(x));
}

/*
 * Převod z formátu plovoucí řádové binární čárky (FP)
 * do formátu pevné řádové binární čárky (FX)
 */
fx fp2fx(double x)
{
    return (fx)(x*(2<<(B-1)));
}

/*
 * Převod z celočíselného formátu (integer)
 * do formátu pevné řádové binární čárky (FX)
 */
fx int2fx(int x)
{
    return (fx)(x<<B);
}

/*
 * Převod z formátu pevné řádové binární čárky (FX)
 * do formátu plovoucí řádové binární čárky (FP)
 */
double fx2fp(fx x)
{
    return (double)x/(2<<(B-1));
}

/*
 * Součet dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_add(fx x, fx y)
{
    return x+y;
}

/*
 * Rozdíl dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_sub(fx x, fx y)
{
    return x-y;
}

/*
 * Součin dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_mul(fx x, fx y)
{
    fx result=(x>>(B/2))*(y>>(B/2));
    return result;
}

/*
 * Podíl dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_div(fx x, fx y)
{
    fx result=x/(y>>(B/2));
    return result<<(B/2);
} 

3. Funkce určená pro naplnění tabulky atans[]

Podobně jako při implementaci algoritmu CORDIC pro hodnoty ve formátu FP, i v případě práce s formátem FX je nutné některé hodnoty předpočítat. Opět se jedná o tabulku atans[] obsahující arkustangenty vybraných úhlů. Kromě toho naplníme i tabulku pows[], ve které jsou uloženy záporné mocniny čísla 2. Jak však uvidíme v dalším textu, není v případě FX výpočtů tato tabulka využita a algoritmus CORDIC se poněkud zjednoduší. Funkce pro výpočet obou zmiňovaných tabulek vypadá následovně:

/*
 * Vytvoření tabulky pro výpočet goniometrických
 * funkcí pomocí algoritmu CORDIC
 */
void fx_create_tables(void)
{
    int i;
    for (i=0; i<MAX_ITER; i++) {
        double p=pow(2.0, -i);
        atans[i]=fp2fx(atan(p));
        pows[i]=fp2fx(p);
    }
} 

To, zda je tabulka korektně naplněna, je možné otestovat velmi jednoduše:

    /* kontrolní výpis tabulky atans[] */
    for (i=0; i<MAX_ITER; i++)
        printf("%d\t%f\n", i, fx2fp(rad2deg(atans[i]))); 

Výsledkem běhu předchozího testu je tabulka hodnot arkustangent úhlů (zhruba) odpovídajících FP verzi:

S takto připravenou tabulkou je možné pokračovat v dalších výpočtech.

4. Výpočet tangenty pomocí algoritmu CORDIC – neoptimalizova­ná verze

Pravděpodobně nejjednodušším výpočtem, který je možné pomocí algoritmu CORDIC provést, je vyjádření tangenty zadaného úhlu. Nejprve si uvedeme neoptimalizovanou verzi, která vznikla v podstatě přímým převodem dříve uvedené implementace určené pro FP reprezentaci. Převod spočívá v náhradě aritmetických funkcí jejich FX ekvivalenty:

/* výpočet funkce tan() pro zadaný úhel delta */
/* (neoptimalizovaná verze) */
fx fx_tan_cordic(fx delta)
{
    int i;
    /* nastavení počátečních podmínek */
    fx x0=fp2fx(1.0);
    fx y0=fp2fx(0.0);
    fx xn;
    if (delta==0) return 0;                 /* ošetření nulového úhlu */
    for (i=0; i<MAX_ITER; i++) {            /* iterační smyčka */
        if (delta<0) {                      /* úhel je záporný => rotace doleva */
            xn=fx_add(x0, fx_mul(y0, pows[i]));
            y0=fx_sub(y0, fx_mul(x0, pows[i]));
            delta=fx_add(delta, atans[i]);
        }
        else {                              /* úhel je kladný => rotace doprava */
            xn=fx_sub(x0, fx_mul(y0, pows[i]));
            y0=fx_add(y0, fx_mul(x0, pows[i]));
            delta=fx_sub(delta, atans[i]);
        }
        x0=xn;
/*        printf("%i\t%+f\t%+f\t%+f\n", i, fx2fp(x0), fx2fp(y0), fx2fp(delta)); */
    }
    if (x0==0)                              /* ošetření tangenty pravého úhlu */
        if (y0<0) return 0;
        else return 0;
    else
        return fx_div(y0,x0);               /* vrátit výsledek operace */
} 

Tuto funkci si můžeme jednoduchým způsobem otestovat, postačí použít smyčku, ve které se počítají tangenty úhlů v rozsahu 0°..89° spolu s vyjádřením absolutních a relativních chyb vzniklých použitím algoritmu CORDIC a bitově omezené FP reprezentace:

    /* výpis tabulky tangent úhlů v rozsahu 0..89° */
    for (i=0; i<90; i++) {                  /* výpočetní smyčka */
        delta=deg2rad(i);                   /* převod úhlu na radiány */
        tanfx=fx_tan_cordic(fp2fx(delta));  /* aplikace algoritmu CORDIC */
        tanval=fx2fp(tanfx);                /* výpočet funkce tan */
        tanerr=fabs(tanval-tan(delta));     /* výpočet absolutních chyb */
                                            /* tisk výsledků */
        printf("%02d\t%14.10f\t%14.10f\t%12.10f\t%8.3f%%\n",
                i,
                tanval,
                tan(delta),
                tanerr,
                tanerr==0 ? 0:100.0*tanerr/tan(delta));
    } 

V tabulce s výsledky si všimněte, že pro některé úhly je vzniklá chyba rovna až deseti procentům. Je to způsobeno nízkým rozlišením použitého FX formátu, zejména hodnot arctat() těch nejmenších úhlů, a chybu již není možné snížit ani zvýšením počtu iterací (ten je nastaven na šestnáct, to odpovídá počtu bitů za binární řádovou čárkou):

5. Výpočet tangenty pomocí algoritmu CORDIC – optimalizovaná ver­ze

Pokud se nad výše uvedenou funkcí fx_tan_cordic() zamyslíme, pravděpodobně zjistíme, že se v ní zbytečně provádí některé aritmetické operace. Zejména se jedná o průběžné násobení koeficienty uloženými v tabulce pows[]. Tato tabulka měla svůj význam při práci s FP formátem, u formátu FX však ztrácí prakticky všechny své výhody, protože násobení zápornou mocninou čísla 2 je vlastně totožné s bitovým posunem doprava. Optimalizovaná funkce pro výpočet tangenty úhlu algoritmem CORDIC vypadá následovně:

/* výpočet funkce tan() pro zadaný úhel delta */
/* (optimalizovaná verze) */
fx fx_tan_cordic_optim(fx delta)
{
    int i;
    /* nastavení počátečních podmínek */
    fx x0=int2fx(1);
    fx y0=0;
    fx xn;
    if (delta==0) return 0;                 /* ošetření nulového úhlu */
    for (i=0; i<MAX_ITER; i++) {            /* iterační smyčka */
        if (delta<0) {                      /* úhel je záporný => rotace doleva */
            xn=fx_add(x0, y0>>i);           /* místo násobení bitový posuv */
            y0=fx_sub(y0, x0>>i);
            delta=fx_add(delta, atans[i]);
        }
        else {                              /* úhel je kladný => rotace doprava */
            xn=fx_sub(x0, y0>>i);
            y0=fx_add(y0, x0>>i);
            delta=fx_sub(delta, atans[i]);
        }
        x0=xn;
    }
    if (x0==0)                              /* ošetření tangenty pravého úhlu */
        if (y0<0) return 0;
        else return 0;
    else
        return fx_div(y0,x0);               /* vrátit výsledek operace */
} 

(Pokud by někdo pocítil potřebu provádět další optimalizace, nabízí se zde použití MMX instrukcí s paralelním výpočtem obou alternativních větví). Za povšimnutí stojí fakt, že se v iterační smyčce používají pouze aritmetické operace součtu, rozdílu a bitového posuvu. Kromě toho se i při inicializaci hodnot x0 a y0 používá pouze celočíselná aritmetika, tj. ve výpočtech se vůbec nevyskytuje volání instrukcí matematického koprocesoru. Tato absence složitých a implementačně náročných operací představuje další z důvodů velké oblíbenosti algoritmu CORDIC v komunitě vývojářů pro malé mikrořadiče a mikroprocesory.

Pokud výše uvedenou funkci fx_tan_cordic_op­tim() podrobíme stejnému testu, jako funkci fx_tan_cordic(), dostaneme následující (do značné míry podobné) výsledky:

6. Výpis celého demonstračního příkladu

V této kapitole je uveden výpis celého demonstračního příkladu, který byl odladěn pro třicetidvoubitové překladače jazyka C. Otestování bylo provedeno pomocí GCC 3.4.2 (uznávám, je to poněkud starší, zato však stabilní verze) a Borland C++ Compileru 5.5 (pro Windows) běžícím v režimu překladače C:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/* počet míst před a za binární řádovou tečkou */
#define A 16
#define B 16

/* Ludolfovo číslo */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* maximální počet iterací při běhu algoritmu */
#define MAX_ITER 16

/* "zesílení" při rotacích */
#define K 0.6073

/* převody mezi stupni a radiány */
#define rad2deg(rad) ((rad)*180.0/M_PI)
#define deg2rad(deg) ((deg)/180.0*M_PI)

/* datový typ, se kterým budeme pracovat */
typedef signed int fx;

/* hlavičky použitých funkcí */
void   fx_print(fx x);
fx     fp2fx(double x);
double fx2fp(fx x);

/* tabulka arkustangentu úhlů */
fx atans[MAX_ITER];

/* tabulka záporných celočíselných mocnin hodnoty 2 */
fx pows[MAX_ITER];

/*
 * Tisk numerické hodnoty uložené ve formátu pevné
 * řádové binární čárky (FX)
 */
void fx_print(fx x)
{
    int i;
    int val=x;                         /* pomocná proměnná pro převod do dvojkové soustavy */
    printf("bin: ");
    for (i=0; i<A+B; i++) {            /* převod na řetězec bitů (do dvojkové soustavy) */
        putchar(!!(val & (1<<(A+B-1)))+'0'); /* výpis hodnoty aktuálně nejvyššího bitu */
        if (i==B-1) putchar('.');      /* po řádové binární čárce vypsat značku */
        val=val<<1;                    /* posun na další (méně významný) bit */
    }

    printf("   hex: %08x   fp: %+11.5f\n", x, fx2fp(x));
}

/*
 * Převod z formátu plovoucí řádové binární čárky (FP)
 * do formátu pevné řádové binární čárky (FX)
 */
fx fp2fx(double x)
{
    return (fx)(x*(2<<(B-1)));
}

/*
 * Převod z celočíselného formátu (integer)
 * do formátu pevné řádové binární čárky (FX)
 */
fx int2fx(int x)
{
    return (fx)(x<<B);
}

/*
 * Převod z formátu pevné řádové binární čárky (FX)
 * do formátu plovoucí řádové binární čárky (FP)
 */
double fx2fp(fx x)
{
    return (double)x/(2<<(B-1));
}

/*
 * Součet dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_add(fx x, fx y)
{
    return x+y;
}

/*
 * Rozdíl dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_sub(fx x, fx y)
{
    return x-y;
}

/*
 * Součin dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_mul(fx x, fx y)
{
    fx result=(x>>(B/2))*(y>>(B/2));
    return result;
}

/*
 * Podíl dvou hodnot uložených ve shodném formátu
 * pevné binární řádové čárky (FX)
 */
fx fx_div(fx x, fx y)
{
    fx result=x/(y>>(B/2));
    return result<<(B/2);
}

/*
 * Vytvoření tabulky pro výpočet goniometrických
 * funkcí pomocí algoritmu CORDIC
 */
void fx_create_tables(void)
{
    int i;
    for (i=0; i<MAX_ITER; i++) {
        double p=pow(2.0, -i);
        atans[i]=fp2fx(atan(p));
        pows[i]=fp2fx(p);
    }
}

/* výpočet funkce tan() pro zadaný úhel delta */
/* (neoptimalizovaná verze) */
fx fx_tan_cordic(fx delta)
{
    int i;
    /* nastavení počátečních podmínek */
    fx x0=fp2fx(1.0);
    fx y0=fp2fx(0.0);
    fx xn;
    if (delta==0) return 0;                 /* ošetření nulového úhlu */
    for (i=0; i<MAX_ITER; i++) {            /* iterační smyčka */
        if (delta<0) {                      /* úhel je záporný => rotace doleva */
            xn=fx_add(x0, fx_mul(y0, pows[i]));
            y0=fx_sub(y0, fx_mul(x0, pows[i]));
            delta=fx_add(delta, atans[i]);
        }
        else {                              /* úhel je kladný => rotace doprava */
            xn=fx_sub(x0, fx_mul(y0, pows[i]));
            y0=fx_add(y0, fx_mul(x0, pows[i]));
            delta=fx_sub(delta, atans[i]);
        }
        x0=xn;
/*        printf("%i\t%+f\t%+f\t%+f\n", i, fx2fp(x0), fx2fp(y0), fx2fp(delta)); */
    }
    if (x0==0)                              /* ošetření tangenty pravého úhlu */
        if (y0<0) return 0;
        else return 0;
    else
        return fx_div(y0,x0);               /* vrátit výsledek operace */
}

/* výpočet funkce tan() pro zadaný úhel delta */
/* (optimalizovaná verze) */
fx fx_tan_cordic_optim(fx delta)
{
    int i;
    /* nastavení počátečních podmínek */
    fx x0=int2fx(1);
    fx y0=0;
    fx xn;
    if (delta==0) return 0;                 /* ošetření nulového úhlu */
    for (i=0; i<MAX_ITER; i++) {            /* iterační smyčka */
        if (delta<0) {                      /* úhel je záporný => rotace doleva */
            xn=fx_add(x0, y0>>i);           /* místo násobení bitový posuv */
            y0=fx_sub(y0, x0>>i);
            delta=fx_add(delta, atans[i]);
        }
        else {                              /* úhel je kladný => rotace doprava */
            xn=fx_sub(x0, y0>>i);
            y0=fx_add(y0, x0>>i);
            delta=fx_sub(delta, atans[i]);
        }
        x0=xn;
    }
    if (x0==0)                              /* ošetření tangenty pravého úhlu */
        if (y0<0) return 0;
        else return 0;
    else
        return fx_div(y0,x0);               /* vrátit výsledek operace */
}


int main(void)
{
    int i;
    fx     tanfx;
    double delta;                       /* úhel, ze kterého se počítá tan */
    double tanval;                      /* vypočtené hodnoty */
    double tanerr;                      /* absolutní chyby */

    fx_create_tables();

    /* kontrolní výpis tabulky atans[] */
    puts("\nHodnoty uložené v tabulce atans[]:");
    for (i=0; i<MAX_ITER; i++)
        printf("%d\t%f\n", i, fx2fp(rad2deg(atans[i])));
        /*printf("%f\n", fx2fp(pows[i])); */

    /* výpočet tan() neooptimalizovanou metodou CORDIC */
    puts("\nVýpočet tan() neoptimalizovanou metodou CORDIC");
    for (i=0; i<90; i++) {                  /* výpočetní smyčka */
        delta=deg2rad(i);                   /* převod úhlu na radiány */
        tanfx=fx_tan_cordic(fp2fx(delta));  /* aplikace algoritmu CORDIC */
        tanval=fx2fp(tanfx);                /* výpočet funkce tan */
        tanerr=fabs(tanval-tan(delta));     /* výpočet absolutních chyb */
                                            /* tisk výsledků */
        printf("%02d\t%14.10f\t%14.10f\t%12.10f\t%8.3f%%\n",
                i,
                tanval,
                tan(delta),
                tanerr,
                tanerr==0 ? 0:100.0*tanerr/tan(delta));
    }

    /* výpočet tan() optimalizovanou metodou CORDIC */
    puts("\nVýpočet tan() optimalizovanou metodou CORDIC");
    for (i=0; i<90; i++) {                  /* výpočetní smyčka */
        delta=deg2rad(i);                   /* převod úhlu na radiány */
        tanfx=fx_tan_cordic_optim(fp2fx(delta));  /* aplikace algoritmu CORDIC */
        tanval=fx2fp(tanfx);                /* výpočet funkce tan */
        tanerr=fabs(tanval-tan(delta));     /* výpočet absolutních chyb */
                                            /* tisk výsledků */
        printf("%02d\t%14.10f\t%14.10f\t%12.10f\t%8.3f%%\n",
                i,
                tanval,
                tan(delta),
                tanerr,
                tanerr==0 ? 0:100.0*tanerr/tan(delta));
    }
    return 0;
} 

7. Obsah poslední části tohoto seriálu

Seriál o numerických formátech s pevnou a plovoucí binární řádovou čárkou se blíží ke svému závěru. V následujícím a současně i posledním dílu si ukážeme implementaci dalších goniometrických funkcí algoritmem CORDIC pracujícím nad FX formátem.