arctan (CORDIC,ベクタモード)

[] [] []


前項では、角度0のベクトル(1,0)をθiずつ回転していって(cosθ,sinθ)を求めたが、逆にベクトル(x,y)が角度0になるように−Δθiずつ回転していくと角度arctan(y/x)が計算できる

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

ついでにrからベクトル(x,y)の長さも求まるが、角度だけ計算するなら不要である。

定数は予め計算しておけばよい。

Δθiの符号は、yがゼロに近づくように決めればよい。すなわちyが正だったらΔθiは負に、yが負になったらΔθiは正にする。

結局、arctan(y/x)=θNを計算する漸化式は次のようになる

角度θの範囲は-0.5π≦θ≦0.5πなので、x0でなければならない(x=0も計算できる)

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

◇◇

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

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

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

◇◇◇

(x,y)=(1,1)の計算例。θ=π/4=0.785398163...

i

xi

yi

θi

atan(2**(-i))

0

1.000000000

1.000000000

0.000000000

0.785398163

1

2.000000000

0.000000000

0.785398163

0.463647609

2

2.000000000

-1.000000000

1.249045772

0.244978663

3

2.250000000

-0.500000000

1.004067109

0.124354995

4

2.312500000

-0.218750000

0.879712115

0.062418810

5

2.326171875

-0.074218750

0.817293305

0.031239833

6

2.328491211

-0.001525879

0.786053471

0.015623729

7

2.328515053

0.034856796

0.770429743

0.007812341

8

2.328787372

0.016665272

0.778242084

0.003906230

9

2.328852470

0.007568447

0.782148314

0.001953123

10

2.328867252

0.003019907

0.784101436

0.000976562

11

2.328870201

0.000745622

0.785077999

0.000488281

12

2.328870566

-0.000391521

0.785566280

0.000244141

13

2.328870661

0.000177051

0.785322139

0.000122070

14

2.328870683

-0.000107235

0.785444209

0.000061035

15

2.328870689

0.000034908

0.785383174

0.000030518

16

2.328870690

-0.000036164

0.785413692

0.000015259

17

2.328870691

-0.000000628

0.785398433

0.000007629

18

2.328870691

0.000017140

0.785390804

0.000003815

19

2.328870691

0.000008256

0.785394618

0.000001907

20

2.328870691

0.000003814

0.785396526

0.000000954

21

2.328870691

0.000001593

0.785397479

0.000000477

22

2.328870691

0.000000482

0.785397956

0.000000238

23

2.328870691

-0.000000073

0.785398195

0.000000119

24

2.328870691

0.000000205

0.785398075

0.000000060

25

2.328870691

0.000000066

0.785398135

0.000000030

26

2.328870691

-0.000000003

0.785398165

0.000000015

27

2.328870691

0.000000031

0.785398150

0.000000007

28

2.328870691

0.000000014

0.785398157

0.000000004

29

2.328870691

0.000000005

0.785398161

0.000000002

30

2.328870691

0.000000001

0.785398163

0.000000001

◇◇◇◇

このアルゴリズムをVivadoHLS合成可能なCコードで書くと次のようになる。

入力c,sは簡単のため1以下としたが、任意の値を入力できるようにap_fixed<W,I>に置き換えてもよい。その際、中間変数x,y,x_tmp,y_tmpはオーバーフローしないようにap_fixed<W+1,I+1>にする。

#include <ap_int.h>

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

 * arctan(s/c)

 *   input range: 0<=c<=1,-1<=s<=1

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

template <int N,int F>

void cordic_vec

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

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

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

    ) {

    /* 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 - vector mode

     * x[0] = c

     * y[0] = s

     * t[0] = 0

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

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

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

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

     */

    ap_fixed<2+F,2> x_tmp;

    ap_fixed<2+F,2> y_tmp;

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

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

    ap_fixed<3+F,3> t = 0;

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

        x_tmp = (x >> i);

        y_tmp = (y >> i);

        if(y>=0) {

            x += y_tmp;

            y -= x_tmp;

            t += atn[i];

        } else {

            x -= y_tmp;

            y += x_tmp;

            t -= atn[i];

        }

    }

    *a = t;

}

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

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> c;  // (i) cos

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

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

    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;

        c = cos(a_val);

        s = sin(a_val);

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

        cout << "cos=" << setw(10) << c << " sin=" << setw(10) << s;

        e = a.to_double() - a_val;

        cout << " a=" << setw(12) << a << " e=" << setw(13) << e;

        cout << endl;

    }

    return 0;

}


[] [] []