sin,cos (CORDIC,回転モード)
CORDICによるsin,cosの計算はsin,cos (CORDIC,回転モード)に書いたとおりで、再掲しておく。
c = cos(πa) s = sin(πa) 角度a:符号つき、整数部1ビット、小数部FAビット (-1≦ a<1) 整数部が1を超える角度は、上位ビット(2の倍数)を切り捨てればよい 余弦cと正弦s:符号つき、整数部2ビット、小数部Fビット (-1≦ s,c ≦1) 小数部はF=FA+2ぐらいがよい。 (∵sinの変化率が最も大きいところでsin(πa)≒πaだからπ倍になる) 誤差:CORDICは収束計算なので、誤差は実際に計算して要確認。 |
VitisHLSで合成可能なCコードは
#include <ap_int.h> #include <cmath> // atan() /*********************************************************** * c = cos(PI*a) s = sin(PI*a) * param range: F=FA+2 or so. * input range: -1<=a<+1 **********************************************************/ template <int F,int FA> void cordic_rot (ap_fixed<1+FA,1> a // (i) angle[radian/PI] ,ap_fixed<2+F,2> *c // (o) cos(PI*a) ,ap_fixed<2+F,2> *s // (o) sin(PI*a) ) { /* atan table */ ap_fixed<1+FA,1,AP_RND> atn[FA+1]; for(int i=-2; i<=FA-2; i++) atn[i+2] = atan(pow(2.0,-i))/3.14159265358979; /* factor 1/K */ 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))/PI * 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<1+FA,1> t = a; for(int i=-2; i<=FA-2; i++) { x_tmp = (x >> i); y_tmp = (y >> i); if(t>=0) { 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 FA=15; // const int F=FA+2; // F=FA+2 or so. const double PI = 3.14159265358979; double a_val,c_val,s_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<1+FA,1> a; // (i) angle[radian/PI] 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を1/M刻みで振ってcos,sinの最大誤差を表示する */ for(int i=-M; i<=M-1; i++) { /* 入力値 */ if(i==-M) { a_val = -1.0; } else { a_val = ((double) i) / ((double) M); } a = a_val; /* 計算実行 */ cordic_rot<F,FA>(a,&c,&s); /* 誤差 */ c_val = c.to_double(); s_val = s.to_double(); ec = c_val - cos(PI*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_val - sin(PI*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); printf("ec_max=%+3.1e at a=%6.3f\n",ec_max,ec_max_a_val); printf("es_min=%+3.1e at a=%6.3f\n",es_min,es_min_a_val); printf("es_max=%+3.1e at a=%6.3f\n",es_max,es_max_a_val); return 0; } |