atan2 (CORDIC,ベクタモード)

[] [] []


前項では、偏角0のベクタ(1,0)θiずつ回転していって(cosθ,sinθ)を求めたが、逆にベクタ(c,s)(r,0)に重なるように−Δθiずつ回転していくと偏角atan2(s,c) が計算できる

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

Δθiの符号は、ベクタ(x,y)(1,0)に近づくようにしたいので、yが正だったらΔθiは負に、yが負になったらΔθiは正にする。

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

上記の式ではベクタ長を一定に保つためcosΔθiをかけているが、今回は偏角が計算できればよいのでこの乗算は省略する。

結局atan2(c,s)=θN計算する漸化式は次のようになる。この漸化式ではベクタ長がK倍に大きくなることに注意。

前項で見たように、初期値をi=0から始めると偏角θ-0.5πθ0.5πしかカバーできないので入力範囲はc0に限定される

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

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

◇◇

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

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

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

またarctanの変化率は最も大きいところでarctan(s/c)=aだから、偏角as/cは同程度のビット幅が必要

すなわち、ベクタ長が短いところではs/cの精度が悪いのでaの精度も悪くなってしまうことに注意

◇◇◇

(c,s)=(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コードで書くと、-0.5πatan2(s,c)0.5πの場合は次のようになる。もし入力c,sの範囲をベクタ長1以下に限定すると、中間変数x,y,x_tmp,y_tmpのベクタ長は入力のK=1/0.6=1.7に収まるので整数部は2ビット幅(-2以上+2未満)で十分である。しかし入力範囲をc,s独立にするならば、入力のベクタ長はc=1,s=1のとき√2になるので、中間変数x,y,x_tmp,y_tmpのベクタ長は√2 K=√2/0.6=2.3、すなわち整数部を3ビット幅(-4以上+4未満)に広げる必要がある。

#include <ap_int.h>

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

* atan2(s,c)

 *  input range: 0<=c<=2**(I-2),-2**(I-2)<=s<=2**(I-2),except for c=s=0

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

template <int I,int F,int FA>   // FA<=30. I+F should be comparable to 2+FA.

void cordic_vec

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

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

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

    ) {

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

        };

    /* 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<I+1+F,I+1> x_tmp;

    ap_fixed<I+1+F,I+1> y_tmp;

    ap_fixed<I+1+F,I+1> x = c;

    ap_fixed<I+1+F,I+1> y = s;

    ap_fixed<2+FA,2>    t = 0.0;

    for(int i=0; i<=FA; 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 I = 2;    // integer part width of c,s

    const int F = 28;   // 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 e,e_min=1.0,e_max=-1.0,e_min_c_val,e_min_s_val,e_max_c_val,e_max_s_val;

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

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

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

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

    const int M = 6 * 1000;

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

        a_val = PI / 2.0 * i / M;

        c = pow(2.0,I-2) * sqrt(2) * cos(a_val);

        s = pow(2.0,I-2) * sqrt(2) * sin(a_val);

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

        e = a.to_double()-atan2(s.to_double(),c.to_double());

        if(e < e_min) {e_min = e; e_min_s_val=s.to_double(); e_min_c_val=c.to_double();}

        if(e > e_max) {e_max = e; e_max_s_val=s.to_double(); e_max_c_val=c.to_double();}

    }

    printf("e_min=%3.1e at s=%6.3f c=%6.3f\n",e_min,e_min_s_val,e_min_c_val);

    printf("e_max= %3.1e at s=%6.3f c=%6.3f\n",e_max,e_max_s_val,e_max_c_val);

    return 0;

}

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

-- I=2 F=28 FA=30

e_min=-1.4e-008 at s=-1.408 c= 0.130

e_max= 1.6e-008 at s= 0.900 c= 1.091

-- I=2 F=29 FA=30

e_min=-7.6e-009 at s= 0.137 c= 1.408

e_max= 8.5e-009 at s= 1.408 c= 0.137

-- I=2 F=30 FA=30

e_min=-5.6e-009 at s=-0.897 c= 1.093

e_max= 6.1e-009 at s= 0.897 c= 1.093

-- I=2 F=31 FA=30

e_min=-5.0e-009 at s=-0.880 c= 1.107

e_max= 4.7e-009 at s= 0.947 c= 1.050

 

◇◇◇◇◇

atan2(s,c)πの場合は次のようになる。入力範囲をc,s独立にするならば、入力のベクタ長はc=1,s=1のとき√2になるので、中間変数x,y,x_tmp,y_tmpのベクタ長は√2 K=√2/0.06=23、すなわち整数部を2ビット幅(-2以上+2未満)から6ビット幅(-32以上+32未満)に広げる必要がある。

#include <ap_int.h>

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

* atan2(s,c)

 *  input range: -2**(I-2)<=c<=2**(I-2),-2**(I-2)<=s<=2**(I-2),except for c=s=0

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

template <int I,int F,int FA>   // FA<=30. I+F should be comparable to FA.

void cordic_vec

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

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

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

    ) {

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

        };

    /* 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<I+4+F,I+4> x_tmp;

    ap_fixed<I+4+F,I+4> y_tmp;

    ap_fixed<I+4+F,I+4> x = c;

    ap_fixed<I+4+F,I+4> y = s;

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

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

        x_tmp = (x >> i);

        y_tmp = (y >> i);

        if(y>=0) {

            x += y_tmp;

            y -= x_tmp;

            t += atn[i+2];

        } else {

            x -= y_tmp;

            y += x_tmp;

            t -= atn[i+2];

        }

    }

    *a = t;

}

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

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

    const int I = 2;    // integer part width of c,s

    const int F = 28;   // 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 e,e_min=1.0,e_max=-1.0,e_min_c_val,e_min_s_val,e_max_c_val,e_max_s_val;

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

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

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

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

    const int M = 12 * 1000;

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

        a_val = PI * i / M;

        c = pow(2.0,I-2) * sqrt(2) * cos(a_val);

        s = pow(2.0,I-2) * sqrt(2) * sin(a_val);

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

        e = a.to_double()-atan2(s.to_double(),c.to_double());

        if(e < e_min) {e_min = e; e_min_s_val=s.to_double(); e_min_c_val=c.to_double();}

        if(e > e_max) {e_max = e; e_max_s_val=s.to_double(); e_max_c_val=c.to_double();}

    }

    printf("e_min=%3.1e at s=%6.3f c=%6.3f\n",e_min,e_min_s_val,e_min_c_val);

    printf("e_max= %3.1e at s=%6.3f c=%6.3f\n",e_max,e_max_s_val,e_max_c_val);

    return 0;

}

誤差を見てみると、ベクタのビット幅I+Fと偏角のビット幅2+FAI+F=FAぐらいがよさそうである。(上記のi=0から始める場合はI+F=2+FAがちょうどよかったが、i=-2から始める場合は2ビット左シフトしているのでFは下位2ビット減らしてもかまわない)

-- I=2 F=26 FA=30

e_min=-9.0e-009 at s=-1.181 c= 0.779

e_max= 1.1e-008 at s= 1.181 c= 0.779

-- I=2 F=27 FA=30

e_min=-6.6e-009 at s=-0.553 c=-1.301

e_max= 7.4e-009 at s= 0.553 c=-1.301

-- I=2 F=28 FA=30

e_min=-5.2e-009 at s=-0.553 c=-1.301

e_max= 5.2e-009 at s= 0.054 c=-1.413

-- I=2 F=29 FA=30

e_min=-4.7e-009 at s=-0.054 c=-1.413

e_max= 5.1e-009 at s= 0.054 c=-1.413


[] [] []