sin,cos (CORDIC,回転モード)

[] [] []


CORDICによるsin,cosの計算はsin,cos (CORDIC,回転モード)に書いたとおりで、再掲しておく。

c = cos(πa)  s = sin(πa)

角度a:符号つき、整数部1ビット、小数部FAビット (-1 a1)

    整数部が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;

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

}


[] [] []