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にすればよい

また、角度の単位をラジアンでなく、12πを1.0として2のべき乗等分にしたい場合は(こうすると2π以上の角度は上位ビットを切り捨てるだけでよく便利)arctanの表に係数をかけておけばよい。

◇◇

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

なので、繰り返し回数Nは、小数部のビット幅以下でよい

実用上はビット幅と繰り返し回数Nを振ってみて決めるとよいかもしれない。

◇◇◇

θ=π/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コードで書くと次のようになる。

#include <ap_int.h>

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

 * cos(a) and sin(a)

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

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

template <int N,int F>

void cordic_rot

    (ap_fixed<2+F,2> a  // (i) angle[radian]

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

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

    ) {

    /* atan table */

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

        {0.785398163,0.463647609,0.244978663,0.124354995

        ,0.062418810,0.031239833,0.015623729,0.007812341

        ,0.003906230,0.001953123,0.000976562,0.000488281

        ,0.000244141,0.000122070,0.000061035,0.000030518

        ,0.000015259,0.000007629,0.000003815,0.000001907

        ,0.000000954,0.000000477,0.000000238,0.000000119

        ,0.000000060,0.000000030,0.000000015,0.000000007

        ,0.000000004,0.000000002,0.000000001

        };

    /* cordic - rotation mode

     * x[0] = 0.607252935

     * 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 = 0.607252935;

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

    ap_fixed<2+F,2> t = a;

    for(int i=0; i<N; 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 N = 29;   // number of iteration (max 31)

    const int F = 30;   // fractional part width

    const double PI = 3.141592653;

    double a_val;

    double e;

    ap_fixed<2+F,2> a;  // (i) angle[radian]

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

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

    cout << endl << "-- N=" << N << " F=" << F << endl;

    const int M = 6;

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

        a_val = PI / 2.0 * i / M;

        a = a_val;

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

        cout << "a=" << setw(9) << a;

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

        cout << " cos=" << setw(12) << c << " e=" << setw(13) << e;

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

        cout << " sin=" << setw(12) << s << " e=" << setw(13) << e;

        cout << endl;

    }

    return 0;

}


[] [] []