atan2 (CORDIC,ベクタモード)

[←]　[↑]　[→]

---

CORDICによるアルゴリズムatan2の計算はatan2 (CORDIC,ベクタモード)に書いたとおりで、再掲しておく。

a = atan(s/c)/π　r = sqrt(s*s+c*c)

ベクタ(c,s)：符号つき、整数部Iビット、小数部Fビット ((c,s)≠(0,0)、sqrt(s*s+c*c)＜2^I-1)

CORDICでは、s,cは独立(正方形)ではなく半径一定(円)で計算するので。

半径r：符号なし、整数部Iビット、小数部Fビット (0≦r＜2^I-1)

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

小数部は(I+F)=(1+FA)-2ぐらいがよい。

(∵atanの変化率が最も大きいところでatan(s/c)/π≒(s/c)/πだから1/π倍になる)

誤差：CORDICは収束計算なので、誤差は実際に計算して要確認。

VitisHLSで合成可能なCコードは

#include <ap_int.h>

#include <cmath>  // atan()

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

* a = atan(s/c)/PI  r = sqrt(s*s+c*c)

 *  param range: (I+F)=(1+FA)-2 or so.

 *  input range: (c,s)!=(0,0),sqrt(s*s+c*c)<2**(I-1)

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

template <int I,int F,int FA>

void cordic_vec

  (ap_fixed<I+F,I> c    // (i) r*cos(PI*a)

  ,ap_fixed<I+F,I> s    // (i) r*sin(PI*a)

  ,ap_fixed<I+F,I> *r   // (o) radius

  ,ap_fixed<1+FA,1> *a  // (o) angle[radian/PI]

  ) {

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

   * 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<1+FA,1>    t = 0.0;

  for(int i=-2; i<=FA-2; 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];

      }

  }

  *r = x * x0;

  *a = t;

}

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

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

  const int I=2,F=9;    //

  const int FA=(I+F)+1;   // (I+F)=(1+FA)-2 or so.

  const double PI = 3.14159265358979;

  double a_val,c_val,s_val;

  double a_expected;

  double er,er_min=1.0,er_max=-1.0

,er_min_c_val,er_min_s_val,er_max_c_val,er_max_s_val;

  double ea,ea_min=1.0,ea_max=-1.0

,ea_min_c_val,ea_min_s_val,ea_max_c_val,ea_max_s_val;

  ap_fixed<I+F,I> c;    // (i) r*cos(PI*a)

  ap_fixed<I+F,I> s;    // (i) r*sin(PI*a)

  ap_fixed<I+F,I> r;    // (o) radius

  ap_fixed<1+FA,1> a;   // (o) angle[radian/PI]

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

  const int M = 12 * 1000;

/* 角度を1/M刻みで振って最大誤差を表示する */

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

/* 入力値 */

    if(i==-M) {

      a_val = -1.0;

      c = pow(2.0,I-1) / sqrt(2) * (-1.0);

      s = pow(2.0,I-1) / sqrt(2) * ( 0.0);

    } else {

      a_val = ((double) i) / ((double) M);

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

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

    }

/* 計算実行 */

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

/* 誤差 */

    c_val = c.to_double();

    s_val = s.to_double();

    er = r.to_double() - sqrt(s_val*s_val+c_val*c_val);

    if(er < er_min)

 {er_min = er; er_min_s_val=s_val; er_min_c_val=c_val;}

    if(er > er_max)

 {er_max = er; er_max_s_val=s_val; er_max_c_val=c_val;}

    a_expected = atan2(s_val,c_val)/PI;

    ea = a.to_double() - a_expected;

// 2.0付近の誤差は、a=1とa=-1の差なので補正する

    if(ea >  1.9) ea = ea - 2.0;

    if(ea < -1.9) ea = ea + 2.0;

    if(ea < ea_min)

 {ea_min = ea; ea_min_s_val=s_val; ea_min_c_val=c_val;}

    if(ea > ea_max)

 {ea_max = ea; ea_max_s_val=s_val; ea_max_c_val=c_val;}

  }

  printf("er_min=%+3.1e at s=%+6.3f c=%+6.3f\n"

,er_min,er_min_s_val,er_min_c_val);

  printf("er_max=%+3.1e at s=%+6.3f c=%+6.3f\n"

,er_max,er_max_s_val,er_max_c_val);

  printf("ea_min=%+3.1e at s=%+6.3f c=%+6.3f\n"

,ea_min,ea_min_s_val,ea_min_c_val);

  printf("ea_max=%+3.1e at s=%+6.3f c=%+6.3f\n"

,ea_max,ea_max_s_val,ea_max_c_val);

  return 0;

}

---

[←]　[↑]　[→]