sin,cos (CORDIC,回転モード)
固定小数点の三角関数はCORDIC(COrdinate Rotational DIgital Computer)という方法で計算できる。このCORDICのアルゴリズムを合成可能なCコードで書いてみよう。
◇
CORDICの基本的なアイデアは、偏角0のベクタ(1,0)をΔθiずつ回転していって(cosθ,sinθ)を求めるというものである。
ここでtanΔθiが2のべき乗になるように回転角を選べば、加減算とシフトだけで回転を計算できる。
Δθiの符号は、Δθiの和がθに一致するようにしたいので、θからΔθiを引いていってゼロに近づくように決めればよい。
また、定数は予め計算してテーブルに持っておけばよい。
結局、(cosθ,sinθ)=(xN,yN)を計算する漸化式は次のようになる。
◇
偏角θのとりうる範囲は、初期値をi=0から始めると次のとおりであり、少なくとも-0.5π≦θ≦0.5πで使える。
-π≦θ≦πで使うには初期値をi=-2から始めればよい。ただし回路規模は大きくなる。
また、偏角の単位をラジアンでなく、1周2πを1.0として2n等分にしたい場合は、arctan2-iのテーブルを2πで割っておけばよい。(こうすると2π以上の偏角は上位ビットを切り捨てるだけでよく便利)
◇◇
演算精度については、漸化式自体は正確なので、入力の量子化誤差と途中結果の丸めだけに依る。
繰り返し回数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 |
◇◇◇◇
アルゴリズムをVitisHLSで合成可能なCコードで書くと、-0.5π≦θ≦0.5πの場合は次のようになる。
|
#include <ap_int.h> /************************************************************** * c=cos(a) s=sin(a) * param range: (2+F) and (2+FA) should be comparable.FA<=30 * input range: -PI/2<=a<=PI/2 *************************************************************/ template <int F,int FA> 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 ) { /* 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); /* 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<FA; i++) { x_tmp = (x >> i); y_tmp = (y >> i); if(t[1+FA]) { 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; /* aをPI/2M刻みで振って,cos,sinの最大誤差を表示する */ 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; } |
誤差は、偏角のビット幅FAを一定にしてベクタのビット幅Fを増やしていくとF=FA+2ぐらいで頭打ちになる。
|
-- F=26 FA=26 ec_min=-1.6e-007 at a=-0.339 ec_max= 2.2e-007 at a= 0.471 es_min=-2.2e-007 at a=-0.029 es_max= 2.2e-007 at a= 0.029 |
|
-- F=27 FA=26 ec_min=-8.3e-008 at a=-0.365 ec_max= 1.0e-007 at a= 0.435 es_min=-8.9e-008 at a= 0.141 es_max= 9.6e-008 at a=-0.059 |
|
-- F=28 FA=26 ec_min=-5.4e-008 at a=-0.365 ec_max= 6.4e-008 at a= 0.435 es_min=-6.2e-008 at a= 0.059 es_max= 6.2e-008 at a=-0.059 |
|
-- F=29 FA=26 ec_min=-4.4e-008 at a=-0.365 ec_max= 5.6e-008 at a=-0.430 es_min=-5.7e-008 at a= 0.134 es_max= 5.5e-008 at a=-0.134 |
|
-- F=30 FA=26 ec_min=-4.3e-008 at a= 0.366 ec_max= 5.3e-008 at a=-0.430 es_min=-5.5e-008 at a= 0.134 es_max= 5.4e-008 at a=-0.134 |
◇◇◇◇◇
-π≦θ≦πの場合は次のようになる。
|
#include <ap_int.h> /************************************************************* * c=cos(a) s=sin(a) * param range: (2+F) and (3+FA) should be comparable.FA<=30 * 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 ) { /* 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)); /* 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<=FA; i++) { x_tmp = (x >> i); y_tmp = (y >> i); if(t[2+FA]) { 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; /* aをPI/M刻みで振って,cos,sinの最大誤差を表示する */ 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; } |
誤差を見ると、偏角のビット幅FAを一定にしてベクタのビット幅Fを増やしていくとF=FA+1ぐらいで頭打ちになる。
|
-- F=29 FA=30 ec_min=-3.5e-008 at a=-0.025 ec_max=+3.5e-008 at a=-0.896 es_min=-3.3e-008 at a=+0.290 es_max=+3.2e-008 at a=-0.771 |
|
-- F=30 FA=30 ec_min=-1.3e-008 at a=+0.785 ec_max=+1.1e-008 at a=-0.804 es_min=-1.2e-008 at a=-0.100 es_max=+1.2e-008 at a=+0.100 |
|
-- F=31 FA=30 ec_min=-7.8e-009 at a=+0.785 ec_max=+6.7e-009 at a=+0.314 es_min=-7.7e-009 at a=-0.032 es_max=+7.7e-009 at a=+0.032 |
|
-- F=32 FA=30 ec_min=-5.1e-009 at a=+0.535 ec_max=+4.5e-009 at a=+0.288 es_min=-4.9e-009 at a=-0.032 es_max=+4.9e-009 at a=+0.032 |
|
-- F=33 FA=30 ec_min=-4.4e-009 at a=+0.531 ec_max=+3.4e-009 at a=+0.488 es_min=-3.7e-009 at a=-0.012 es_max=+3.7e-009 at a=+0.012 |