Algoritmus CORDIC v FX formátu a goniometrické funkce

Obsah

1. Výpočet goniometrických funkcí sin() a cos()
2. Funkce pro výpočet sin() algoritmem CORDIC
3. Vliv počtu iterací na přesnost výpočtů funkce sin()
4. Funkce pro výpočet cos() algoritmem CORDIC
5. Vliv počtu iterací na přesnost výpočtů funkce cos()
6. Literatura
7. Odkazy na Internetu
8. Slovo závěrem

1. Výpočet goniometrických funkcí sin() a cos()

Jak je již z názvu první kapitoly patrné, dnes se budeme zabývat výpočty goniometrických funkcí pomocí algoritmu CORDIC, přičemž všechny výpočty budou prováděny ve formátu pevné řádové binární čárky. Konkrétně se bude jednat o formát čísel se znaménkem, kde počet bitů před binární řádovou čárkou i za binární řádovou čárkou je shodný a dosahuje hodnoty šestnácti bitů: A=B=16. Na začátku si pro zopakování řekněme, jakým způsobem se algoritmus CORDIC používá pro výpočet funkcí sinus a kosinus. Nejprve je vypočtena tabulka úhlů (tj. arkustangenty záporných druhých mocnin dvojky) naprosto stejným způsobem jako v demonstračním příkladu uvedeném minule. Po vytvoření této tabulky je již možné CORDIC spustit. Počáteční souřadnice vektoru r jsou nastaveny na hodnotu [1, 0] (samozřejmě až po převodu na formát pevné řádové binární čárky, což je ovšem pro hodnoty 1 a 0 triviální). Vektor je posléze v iterační smyčce rotován tak dlouho, dokud neproběhne daný počet iterací.

Úhel vektoru r se přitom z počáteční hodnoty 0 neustále přibližuje k zadanému úhlu δ, jelikož se v iterační smyčce zadaný úhel adaptivně buď zmenšuje či zvětšuje o hodnotu uloženou v tabulce atans[] (zde se jedná o posloupnost, ve které se čísla neustále zmenšují, proto je konvergence ke správnému výsledku zaručena). Výsledek, tj. hodnoty funkcí sinus a kosinus, je uložen v nových souřadnicích vektoru r (v našem případě vynásobený o konstantu K_fx, což je ve skutečnosti konstanta F_float převedená pro snazší výpočty do formátu FX), a to z toho důvodu, že vektor rotoval na jednotkové kružnici a souřadnice jakéhokoli bodu ležícího na jednotkové kružnici přímo odpovídají hodnotám sinu a kosinu úhlu tohoto bodu a počátku souřadné soustavy počítaného od kladné horizontální poloosy. Nejprve si ukažme funkci, která vypočte tabulku úhlů pro záporné druhé mocniny čísla dvě (i tuto funkci jsme si uváděli minule):

#define MAX_ITER dostatečně_velké_číslo

/*
 * 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);
    }
} 

Obsah tabulky si můžeme pro jistotu vypsat na standardní výstup:

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

Dále si uveďme výpočet konstanty K uložené ve formátu pevné řádové binární čárky:

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

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

/* převod K do formátu FX */
static fx K_fx=(fx)(K_float*(2<<(B-1))); 

2. Funkce pro výpočet sin() algoritmem CORDIC

Výpočet goniometrické funkce sin() pomocí algoritmu CORDIC je při použití formátu pevné řádové binární čárky velmi snadný a současně i přímočarý. Dokonce je možné říci, že výpočet je oproti verzi používající plovoucí řádovou čárku jednodušší, a to z toho důvodu, že se místo pomocné tabulky s druhými mocninami čísla dvě přímo používá operace pro aritmetický posuv doprava (ten je podporován přímo centrální procesorovou jednotkou, jedná se tedy o znatelné zrychlení celého výpočtu). Při implementaci algoritmu CORDIC v programovacím jazyku C může funkce pro výpočet sinu vypadat následovně (připomeňme si, že výsledek je uložen v y-ové složce orotovaného vektoru r):

/* výpočet funkce sin() pro zadaný úhel delta */
fx fx_sin_cordic_optim(fx delta)
{
    int i;
    static fx K_fx=(fx)(K_float*(2<<(B-1)));
    /* nastavení počátečních podmínek */
    fx x0=int2fx(1);
    fx y0=0;
    fx xn;
    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;
    }
    return fx_mul(y0, K_fx);                /* opravit "zesílení" výsledku */
} 

Výše uvedenou Céčkovskou funkci fx_sin_cordic_op­tim() si můžeme otestovat na výpočtu úhlů ležících v prvním kvadrantu, tj. od 0° do 90°. Kromě vypočtených hodnot je u každého úhlu uvedena i hodnota vypočtená pomocí FPU, dále pak absolutní a relativní chyba. Tučně jsou zvýrazněny ty relativní chyby, které jsou menší než jedno procento:

Výpočet funkce sin() optimalizovanou metodou CORDIC

Statistické vlastnosti výpočtu

3. Vliv počtu iterací na přesnost výpočtů funkce sin()

V předchozí kapitole jsme si uvedli výsledky výpočtů funkce sinus pro šestnáct iterací. Proč jsem však zvolil zrovna tuto hodnotu? Stanovit optimální počet iterací je poněkud ošemetné. Malý počet iterací má za následek zvýšení chyb při výpočtu, výpočty však budou prováděny značně rychleji. Naopak příliš velký počet iterací může výpočet zpomalit, a to bez vlivu na přesnost výsledku, protože při použití FX aritmetiky je omezený pouze počet bitů za binární řádovou čárkou. Jako pravidlo je možné uvést, že počet iterací by měl být zhruba shodný s počtem významových bitů umístěných za binární řádovou čárkou. Dokonce se o tom můžeme přesvědčit při použití mírně upravené funkce pro výpočet sinu algoritmem CORDIC:

/* výpočet funkce sin() pro zadaný úhel delta */
/* a maximální počet iterací                  */
fx fx_sin_cordic_optim_iter(fx delta, int iter)
{
    int i;
    static fx K_fx=(fx)(K_float*(2<<(B-1)));
    /* nastavení počátečních podmínek */
    fx x0=int2fx(1);
    fx y0=0;
    fx xn;
    for (i=0; i<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;
    }
    return fx_mul(y0, K_fx);                /* opravit "zesílení" výsledku */
} 

Tuto funkci využijeme v následujícím testovacím programu, který vypočte statistické informace o chybách pro maximální počet iterací jdoucí od 1 do 32:

int main(void)
{
    int i;
    fx     sinfx;
    double delta;                           /* úhel, ze kterého se funkce počítá */
    double value;                           /* vypočtené hodnoty */
    double abs_err;                         /* absolutní chyby */
    double rel_err;                         /* relativní chyby */
    double sum_abs_err=0.0;                 /* součet chyb */
    double sum_rel_err=0.0;                 /* součet chyb */
    int iter;                               /* počitadlo iterací */

    fx_create_tables();

    /* výpočet sin() optimalizovanou metodou CORDIC */
    puts("\n<h2>Vliv počtu iterací na výpočet sin() optimalizovanou metodou CORDIC</h2>\n");
    puts("<table>");
    puts("<tr>"\
         "<th><strong>Počet iterací:</strong></th>"\
         "<th><strong>Součet absolutních chyb:</strong></th>"\
         "<th><strong>Součet relativních chyb:</strong></th>"\
         "<th><strong>Průměrná absolutní chyba:</strong></th>"\
         "<th><strong>Průměrná relativní chyba:</strong></th></tr>");
    for (iter=1; iter<=32; iter++) {
        sum_abs_err=0.0;
        sum_rel_err=0.0;
        for (i=0; i<=90; i++) {                 /* výpočetní smyčka */
            delta=deg2rad(i);                   /* převod úhlu na radiány */
            sinfx=fx_sin_cordic_optim_iter(fp2fx(delta), iter);  /* aplikace algoritmu CORDIC */
            value=fx2fp(sinfx);                 /* výpočet funkce sin */
            abs_err=fabs(value-sin(delta));     /* výpočet absolutních chyb */
            rel_err=sin(delta)==0 ? 0:100.0*abs_err/sin(delta);
            sum_abs_err+=abs_err;
            sum_rel_err+=rel_err;
        }
        printf("<tr><td>%d</td><td>%5.3f</td><td>%5.3f%%</td><td>%5.3f</td><td>%5.3f%%</td></tr>\n",
               iter, sum_abs_err, sum_rel_err, sum_abs_err/91.0, sum_rel_err/91.0);
    }
    puts("</table>");
    return 0;
} 

Výsledky běhu předešlého testovacího příkladu jsou velmi zajímavé, zejména při pohledu na zlom okolo hranice čtrnácti iterací:

Vliv počtu iterací na výpočet sin() optimalizovanou metodou CORDIC

4. Funkce pro výpočet cos() algoritmem CORDIC

Goniometrická funkce kosinus se při použití algoritmu CORDIC počítá prakticky stejným způsobem, jako funkce sinus. Jediný rozdíl spočívá v tom, že po dokončení rotace vektoru r je hodnota sinu uložena v jeho y-ové složce, zatímco hodnota kosinu ve složce x-ové. Zdrojový tvar céčkovské funkce na výpočet kosinu je tedy nepatrně odlišný v posledním řádku:

/* výpočet funkce cos() pro zadaný úhel delta */
fx fx_cos_cordic_optim(fx delta)
{
    int i;
    static fx K_fx=(fx)(K_float*(2<<(B-1)));
    /* nastavení počátečních podmínek */
    fx x0=int2fx(1);
    fx y0=0;
    fx xn;
    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;
    }
    return fx_mul(x0, K_fx);                /* opravit "zesílení" výsledku */
} 

Výpočet funkce cos() optimalizovanou metodou CORDIC

Výsledky výpočtu funkce kosinus pro šestnáct iterací algoritmu CORDIC jsou zobrazeny v následující tabulce. Opět jsou zvýrazněny ty úhly, pro něž je relativní chyba menší než jedno procento.

5. Vliv počtu iterací na přesnost výpočtů funkce cos()

Podobně jako u funkce sinus, i u kosinu je možné zjistit vliv maximálního počtu iterací algoritmu CORDIC na přesnost výpočtů funkce. Jelikož jsme použili FX formát s šestnácti bity umístěnými za binární řádovou čárkou, je možné předpokládat, že optimální bude použití šestnácti iterací, o čemž se můžeme přesvědčit po překladu a spuštění následujícího testovacího příkladu.

Vliv počtu iterací na výpočet cos() optimalizovanou metodou CORDIC

int main(void)
{
    int i;
    fx     cosfx;
    double delta;                           /* úhel, ze kterého se funkce počítá */
    double value;                           /* vypočtené hodnoty */
    double abs_err;                         /* absolutní chyby */
    double rel_err;                         /* relativní chyby */
    char  *zvyr1, *zvyr2;                   /* ukazatele na konstantní řetězce pro */
                                            /* generování HTML */

    fx_create_tables();

    puts("\n<h2>Výpočet funkce cos() optimalizovanou metodou CORDIC</h2>\n");
    puts("<table>");
    for (i=0; i<=90; i++) {                 /* výpočetní smyčka */
        delta=deg2rad(i);                   /* převod úhlu na radiány */
        cosfx=fx_cos_cordic_optim(fp2fx(delta));  /* aplikace algoritmu CORDIC */
        value=fx2fp(cosfx);                 /* výpočet funkce cos */
        abs_err=fabs(value-cos(delta));     /* výpočet absolutních chyb */
        rel_err=cos(delta)<=1e-10 ? 0:100.0*abs_err/cos(delta);
        if (rel_err<=1.0) {
            zvyr1="<strong>";
            zvyr2="</strong>";
        }
        else {
            zvyr1="";
            zvyr2="";
        }
        printf("<tr><td>%02d</td><td>%5.3f</td><td>%5.3f%</td><td>%5.3f</td><td>%s%5.3f%%%s</td></tr>\n",
           i, value, cos(delta), abs_err, zvyr1, rel_err, zvyr2);
    }
    puts("</table>");
    return 0;
} 

Výsledky běhu testovacího příkladu odpovídají hodnotám zjištěným u funkce sinus:

6. Literatura

1. Yates Randy: Fixed-Point Arithmetic: An Introduction,
   Digital Sound Labs, March 3, 2001
2. Hook Brian: An Introduction to Fixed Point Math,
   Game Design and Review, 2003
3. Andraka, Ray: „A survey of CORDIC algorithms for FPGA based computers“,
   ACM, 1998
4. Despain, A.M.:„Fourier Transform Computations Using CORDIC Iterations“,
   IEEE Transactions on Computers, Volume 23, strany 993–1001, 1974
5. Mazenc C., Merrheim, X., Muller, J.M.: „Computing Functions Arccos and Arcsin Using CORDIC“,
   IEEE Transactions on Computers, Volume 42, strany 118–122, 1993
6. Volder, Jack: „Binary computation algorithms for coordinate rotation and function generation“,
   Convair Report IAR-1, 1956
7. Volder, Jack: „The CORDIC Trigonometric Computing Technique“,
   IRE Transactions on Electronic Computing, Vol EC-8, strany 330–334, 1959
8. NVIDIA Corporation: „Floating-Point Specials on the GPU“,
   2005

7. Odkazy na Internetu

1. P. Mikulec, M. Vojtíšek: Procesor IBM RS 6000,
   http://petam.chy­trak.cz/skola/RS6000
2. Wikipedia: Fixed-point arithmetic,
   http://en.wiki­pedia.org/wiki/Fi­xed-point_arithmetic
3. Wikipedia: Floating point,
   http://en.wiki­pedia.org/wiki/Flo­ating_point
4. Wikipedia: IEEE floating-point standard,
   http://en.wiki­pedia.org/wiki/I­EEE_Floating_Po­int_Standard
5. Grant R. Griffin: CORDIC FAQ,
   http://www.dspgu­ru.com/info/faq­s/cordic.htm
6. Andraka Consulting Group, Inc.: What is all this CORDIC stuff anyhow?,
   http://www.an­draka.com/cor­dic.htm
7. Cyliax Ingo: CORDIC (COordinate Rotation DIgital Computer), the swiss army knife for computing math functions…
   http://www.ee­.ualberta.ca/cou­rses/ee401/mi­croboard/cordic_CCin­k.html

8. Slovo závěrem

V tomto desetidílném seriálu jsme si přiblížili způsob použití dvou základních formátů podmnožiny racionálních čísel v operační paměti počítače. Jednalo se o široce používaný a podporovaný formát pohyblivé řádové binární čárky a již méně známý, ale stále užitečný a nenahraditelný formát pevné řádové binární čárky. Popsali jsme si vlastnosti obou formátů, včetně předností a záporů při jejich použití v různých typech úloh, jakými jsou počty s měnou, zpracování signálů, výpočty goniometrických funkcí atd. Také jsme si ukázali, jakým způsobem jsou v obou zmiňovaných formátech implementovány základní aritmetické operace (součet, rozdíl, součin, podíl), ale i operace složitější – výpočet druhé odmocniny, vyčíslení goniometrických funkcí, výpočet délky vektoru atd. V neposlední řadě byl v několika demonstračních příkladech prakticky ukázán algoritmus CORDIC, který má v mnoha systémech i výpočetně náročných aplikacích svoji nezastupitelnou­ úlohu.