atan2 (CORDIC,ベクタモード)
前項では、偏角0のベクタ(1,0)をθiずつ回転していって(cosθ,sinθ)を求めたが、逆にベクタ(c,s)が(r,0)に重なるように−Δθiずつ回転していくと偏角atan2(s,c) が計算できる。
前項同様、tanΔθiが2のべき乗になるように回転角を選べば、加減算とシフトだけで回転を計算できる。
Δθiの符号は、ベクタ(x,y)が(1,0)に近づくようにしたいので、yが正だったらΔθiは負に、yが負になったらΔθiは正にする。また、定数は予め計算してテーブルに持っておけばよい。
上記の式ではベクタ長を一定に保つためcosΔθiをかけているが、今回は偏角が計算できればよいのでこの乗算は省略する。
結局、atan2(c,s)=θNを計算する漸化式は次のようになる。この漸化式ではベクタ長がK倍に大きくなることに注意。
◇
前項で見たように、初期値をi=0から始めると偏角θは-0.5π≦θ≦0.5πしかカバーできないので入力範囲はc≧0に限定される。-π≦θ≦πで使うには初期値をi=-2から始めればよい。ただし回路規模は大きくなる。
また、偏角の単位をラジアンでなく、1周2πを1.0として2n等分にしたい場合は、arctan2-iのテーブルを2πで割っておけばよい。(こうすると2π以上の偏角は上位ビットを切り捨てるだけでよく便利)
◇◇
演算精度については、漸化式自体は正確なので、入力の量子化誤差と途中結果の丸めだけに依る。
繰り返し回数Nは、偏角の刻みが
となるから、偏角の小数部ビット幅まで繰り返せば十分である。
またarctanの変化率は最も大きいところでarctan(s/c)=aだから、偏角aとs/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 |
◇◇◇◇
アルゴリズムをVitisHLSで合成可能な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> /************************************************************* * a = atan2(s,c) * param range: I+F and 2+FA should be comparable. FA<=30 * input range: 0<=c<=2**(I-2),-2**(I-2)<=s<=2**(I-2) * (c,s)!=(0,0) ************************************************************/ template <int I,int F,int 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] ) { /* atan table */ ap_fixed<1+FA,1,AP_RND> atn[31]; for(int i=0; i<=30; i++) atn[i] = atan(pow(2.0,-i)); /* 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[I+F]) { 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; } |
誤差を見ると、偏角のビット幅FAを一定してベクタのビット幅I+Fを増やしていくとI+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> /************************************************************** * a = atan2(s,c) * param range: I+F and 3+FA should be comparable. FA<=30 * input range: -2**(I-2)<=c<=2**(I-2),-2**(I-2)<=s<=2**(I-2) * (c,s)!=(0,0) **************************************************************/ template <int I,int F,int 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] ) { /* 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)); /* 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[I+3+F]) { 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; } |
誤差を見ると、偏角のビット幅FAを一定してベクタのビット幅I+Fを増やしていくとI+F=FAぐらいで頭打ちになる。(上記のi=0から始める場合はI+F=2+FAがちょうどよかったが、i=-2から始める場合は2ビット左シフトしているのでFは下位2ビット減らしてもかまわない)
|
-- I=2 F=25 FA=30 e_min=-1.5e-008 at s=+1.323 c=-0.500 e_max=+1.8e-008 at s=+1.112 c=+0.873 |
|
-- I=2 F=26 FA=30 e_min=-8.8e-009 at s=+1.345 c=-0.437 e_max=+1.0e-008 at s=+1.181 c=+0.779 |
|
-- I=2 F=27 FA=30 e_min=-6.2e-009 at s=+1.323 c=-0.500 e_max=+6.4e-009 at s=+0.553 c=-1.301 |
|
-- I=2 F=28 FA=30 e_min=-4.9e-009 at s=+1.345 c=-0.437 e_max=+4.8e-009 at s=-1.209 c=+0.733 |
|
-- I=2 F=29 FA=30 e_min=-4.3e-009 at s=+1.389 c=-0.268 e_max=+4.2e-009 at s=+0.054 c=-1.413 |