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

[] [] []


固定小数点の三角関数はCORDIC(COrdinate Rotational DIgital Computer)という方法で計算できる。このCORDICのアルゴリズムを合成可能なCコードで書いてみよう

CORDICの基本的なアイデアは、偏角0のベクタ(1,0)Δθiずつ回転していって(cosθ,sinθ)を求めるというものである

ここでtanΔθi2のべき乗になるように回転角を選べば、加減算とシフトだけで回転を計算できる。

Δθiの符号は、Δθiの和がθに一致するようにしたいので、θからΔθiを引いていってゼロに近づくように決めればよい。

また、定数は予め計算してテーブルに持っておけばよい。

結局、(cosθ,sinθ)=(xN,yN)を計算する漸化式は次のようになる

偏角θのとりうる範囲は、初期値をi=0から始めると次のとおりであり、少なくとも-0.5πθ0.5πで使える

θπで使うには初期値i=-2から始めればよい。ただし回路規模は大きくなる。

また、偏角の単位をラジアンでなく、11.0として2n等分にしたい場合は(こうすると以上の偏角は上位ビットを切り捨てるだけでよく便利)arctan2-iのテーブルをで割っておけばよい。

◇◇

演算精度については、漸化式自体は正確なので、入力の量子化誤差と途中結果の丸めだけに依る

繰り返し回N、偏角の刻みが

となるから、偏角の小数部ビット幅まで繰り返せば十分である

またsin(a)の変化率は最も大きいところでsin(a)=aだから、偏角aとベクタ(x,y)は同程度のビット幅が必要。実際に値を計算してみて決めるのがよいかもしれない。

◇◇◇

θ=π/4の計算例。cosθ=sinθ=1/√2=0.707106781...

i

xi

yi

θi

atan(2**(-i))

0

0.607252935

0.000000000

0.785398163

0.785398163

1

0.607252935

0.607252935

0.000000000

0.463647609

2

0.303626468

0.910879403

-0.463647609

0.244978663

3

0.531346318

0.834972786

-0.218668946

0.124354995

4

0.635717916

0.768554496

-0.094313951

0.062418810

5

0.683752572

0.728822126

-0.031895141

0.031239833

6

0.706528264

0.707454858

-0.000655308

0.015623729

7

0.717582246

0.696415354

0.014968421

0.007812341

8

0.712141501

0.702021465

0.007156080

0.003906230

9

0.709399230

0.704803268

0.003249850

0.001953123

10

0.708022661

0.706188813

0.001296727

0.000976562

11

0.707333023

0.706880242

0.000320165

0.000488281

12

0.706987867

0.707225619

-0.000168116

0.000244141

13

0.707160529

0.707053015

0.000076024

0.000122070

14

0.707074219

0.707139338

-0.000046046

0.000061035

15

0.707117380

0.707096182

0.000014989

0.000030518

16

0.707095801

0.707117761

-0.000015529

0.000015259

17

0.707106590

0.707106972

-0.000000270

0.000007629

18

0.707111985

0.707101577

0.000007360

0.000003815

19

0.707109288

0.707104274

0.000003545

0.000001907

20

0.707107939

0.707105623

0.000001638

0.000000954

21

0.707107265

0.707106298

0.000000684

0.000000477

22

0.707106928

0.707106635

0.000000207

0.000000238

23

0.707106759

0.707106803

-0.000000031

0.000000119

24

0.707106843

0.707106719

0.000000088

0.000000060

25

0.707106801

0.707106761

0.000000028

0.000000030

26

0.707106780

0.707106782

-0.000000001

0.000000015

27

0.707106791

0.707106772

0.000000013

0.000000007

28

0.707106785

0.707106777

0.000000006

0.000000004

29

0.707106783

0.707106780

0.000000002

0.000000002

30

0.707106781

0.707106781

0.000000000

0.000000001

◇◇◇◇

アルゴリズムVivadoHLSで合成可能なCコードで書くと、-0.5πθ0.5πの場合は次のようになる。

#include <ap_int.h>

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

 * cos(a) and sin(a)

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

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

template <int F,int FA> // FA<=30. F and FA should be comparable.

void cordic_rot

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

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

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

    ) {

    /* atn[i] = atan(pow(2.0,-i)) */

    ap_fixed<1+FA,1,AP_RND> atn[31] =

        {0.7853981634,0.4636476090,0.2449786631,0.1243549945

        ,0.0624188100,0.0312398334,0.0156237286,0.0078123411

        ,0.0039062301,0.0019531225,0.0009765622,0.0004882812

        ,0.0002441406,0.0001220703,0.0000610352,0.0000305176

        ,0.0000152588,0.0000076294,0.0000038147,0.0000019073

        ,0.0000009537,0.0000004768,0.0000002384,0.0000001192

        ,0.0000000596,0.0000000298,0.0000000149,0.0000000075

        ,0.0000000037,0.0000000019,0.0000000009

        };

    /* initial value x0 */

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

    /* 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<2+FA,2> t = a;

    for(int i=0; i<F; i++) {

        x_tmp = (x >> i);

        y_tmp = (y >> i);

        if(t>=0) {

            x -= y_tmp;

            y += x_tmp;

            t -= atn[i];

        } else {

            x += y_tmp;

            y -= x_tmp;

            t += atn[i];

        }

    }

    *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<2+FA,2> 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 = 6 * 1000;

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

        a_val = PI / 2.0 * 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;

}

誤差を見てみると、ベクタの小数部ビット幅Fと偏角の小数部ビット幅FAF=FA+2ぐらいがよさそうである。

-- F=30 FA=30

ec_min=-1.1e-008 at a= 0.326

ec_max= 1.4e-008 at a=-0.203

es_min=-1.1e-008 at a= 0.174

es_max= 1.1e-008 at a= 0.029

-- F=31 FA=30

ec_min=-6.0e-009 at a= 0.365

ec_max= 6.9e-009 at a= 0.490

es_min=-6.2e-009 at a=-0.334

es_max= 6.2e-009 at a= 0.334

-- F=32 FA=30

ec_min=-3.4e-009 at a= 0.365

ec_max= 4.4e-009 at a= 0.471

es_min=-3.9e-009 at a= 0.135

es_max= 3.9e-009 at a=-0.135

-- F=33 FA=30

ec_min=-2.9e-009 at a=-0.287

ec_max= 2.8e-009 at a= 0.471

es_min=-3.2e-009 at a= 0.174

es_max= 3.2e-009 at a=-0.174

 

◇◇◇◇◇

θ0πの場合は次のようになる。

#include <ap_int.h>

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

 * cos(a) and sin(a)

 *   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

    ) {

    /* atn[i+2] = atan(pow(2.0,-i)) */

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

        {1.3258176637,1.1071487178

        ,0.7853981634,0.4636476090,0.2449786631,0.1243549945

        ,0.0624188100,0.0312398334,0.0156237286,0.0078123411

        ,0.0039062301,0.0019531225,0.0009765622,0.0004882812

        ,0.0002441406,0.0001220703,0.0000610352,0.0000305176

        ,0.0000152588,0.0000076294,0.0000038147,0.0000019073

        ,0.0000009537,0.0000004768,0.0000002384,0.0000001192

        ,0.0000000596,0.0000000298,0.0000000149,0.0000000075

        ,0.0000000037,0.0000000019,0.0000000009

        };

    /* 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>=0) {

            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;

    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;

}

誤差を見てみると、ベクタの小数部ビット幅Fと偏角の小数部ビット幅FAF=FA+2ぐらいがよさそうである。

-- F=30 FA=30

ec_min=-1.1e-008 at a= 0.559

ec_max= 1.3e-008 at a=-0.804

es_min=-1.1e-008 at a= 0.025

es_max= 1.1e-008 at a=-0.025

-- F=31 FA=30

ec_min=-6.0e-009 at a= 0.600

ec_max= 8.2e-009 at a= 0.468

es_min=-6.1e-009 at a= 0.290

es_max= 6.1e-009 at a=-0.290

-- F=32 FA=30

ec_min=-4.0e-009 at a= 0.556

ec_max= 5.1e-009 at a= 0.293

es_min=-5.0e-009 at a=-0.978

es_max= 4.8e-009 at a= 0.909

-- F=33 FA=30

ec_min=-3.3e-009 at a=-0.561

ec_max= 4.1e-009 at a=-0.293

es_min=-4.5e-009 at a=-0.889

es_max= 4.5e-009 at a= 0.889


[] [] []