sin,cos (CORDIC,回転モード)

[←]　[↑]　[→]

---

CORDICによるsin,cosの計算はsin,cos (CORDIC,回転モード)に書いたとおりで、再掲しておく。

偏角a：符号つき、整数部3ビット、小数部FAビット

　　　　-π≦ a ≦π

　　　　FAは定数表atn[i]の精度に収める(FA≦30)

　　　　単位をラジアンでなく1周2πを1.0として2ⁿ等分にする場合は

　　　　定数表atn[i]を2πで割っておけばよい

余弦cと正弦s：符号つき、整数部2ビット、小数部Fビット

　　　　-1≦ s,c ≦1

誤差：CORDICは収束計算なので、誤差は実際に計算して要確認。

　　　またsin,cosの変化率は大きいところでaと同じだから、

　　　同じくらいの精度(FA≒F)にする必要がある。

VitisHLSで合成可能なCコードは

#include <ap_int.h>

#include <cmath>  // atan()

/*************************************************************

 * c=cos(a) s=sin(a)

 *   param range: FA<=30. F should be comparable to FA.

 *   input range: -PI<=a<=PI

 ************************************************************/

template <int F,int FA>

void cordic_rot

    (ap_fixed<3+FA,3> a // (i) argument[radian]

    ,ap_fixed<2+F,2> *c // (o) cos

    ,ap_fixed<2+F,2> *s // (o) sin

    ) {

/* atan table */

    ap_fixed<2+FA,2,AP_RND> atn[33];

    for(int i=-2; i<=30; i++) atn[i+2] = atan(pow(2.0,-i));

/* initial value x0 */

    ap_fixed<1+F,1,AP_RND> x0 = 0.0658658286;

/* cordic - rotation mode

     * x[0] = x0

     * y[0] = 0

     * t[0] = a

     * x[i+1] = x[i] - sign(t[i])*y[i]*2**(-i)

     * y[i+1] = y[i] + sign(t[i])*x[i]*2**(-i)

     * t[i+1] = t[i] - sign(t[i])*atan(2**(-i))

     * where sign(t) = (t>=0)? +1 : -1

     */

    ap_fixed<2+F,2>  x_tmp;

    ap_fixed<2+F,2>  y_tmp;

    ap_fixed<2+F,2>  x = x0;

    ap_fixed<2+F,2>  y = 0.0;

    ap_fixed<3+FA,3> t = a;

    for(int i=-2; i<=F; i++) {

        x_tmp = (x >> i);

        y_tmp = (y >> i);

        if(t[2+FA]) {

            x += y_tmp;

            y -= x_tmp;

            t += atn[i+2];

        } else {

            x -= y_tmp;

            y += x_tmp;

            t -= atn[i+2];

        }

    }

    *c = x;

    *s = y;

}

このコードは次のテストベンチで検証してある。

int main(int argc, char *argv[]) {

    const int F = 32;  // fractional part width of c,s

    const int FA = 30; // fractional part width of a (max 30)

    const double PI = 3.14159265358979;

    double a_val;

    double ec,ec_min=1.0,ec_max=-1.0,ec_min_a_val,ec_max_a_val;

    double es,es_min=1.0,es_max=-1.0,es_min_a_val,es_max_a_val;

    ap_fixed<3+FA,3> a; // (i) argument[radian]

    ap_fixed<2+F,2> c;  // (o) cos

    ap_fixed<2+F,2> s;  // (o) sin

    printf("-- F=%-d FA=%-d\n",F,FA);

    const int M = 12 * 1000;

    /* aをPI/M刻みで振ってcos,sinの最大誤差を表示する */

    for(int i=-M; i<=M; i++) {

        a_val = PI * i / M;

        a = a_val;

        cordic_rot<F,FA>(a,&c,&s);

        ec = c.to_double()-cos(a_val);

        if(ec < ec_min) {ec_min = ec; ec_min_a_val = a_val;}

        if(ec > ec_max) {ec_max = ec; ec_max_a_val = a_val;}

        es = s.to_double()-sin(a_val);

        if(es < es_min) {es_min = es; es_min_a_val = a_val;}

        if(es > es_max) {es_max = es; es_max_a_val = a_val;}

    }

    printf("ec_min=%+3.1e at a=%6.3f\n",ec_min,ec_min_a_val/PI);

    printf("ec_max=%+3.1e at a=%6.3f\n",ec_max,ec_max_a_val/PI);

    printf("es_min=%+3.1e at a=%6.3f\n",es_min,es_min_a_val/PI);

    printf("es_max=%+3.1e at a=%6.3f\n",es_max,es_max_a_val/PI);

    return 0;

}

---

[←]　[↑]　[→]