atan2 (CORDIC,ベクタモード)

[←]　[↑]　[→]

---

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

横cと縦s：符号つき、整数部Iビット、小数部Fビット

　　　　　-2^I-2 ≦ s,c ≦2^I-2

偏角a：符号つき、整数部3ビット、小数部FAビット

　　　　-π≦ a ≦π

　　　　FAは定数表atan[i]の精度に収める(FA≦30)

　　　　単位をラジアンでなく1周2πを1.0として2ⁿ等分にする場合は

　　　　定数表atn[i]を2πで割っておけばよい

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

　　　またatanの変化率は大きいところでs/cと同じだから、

　　　同じくらいの精度(FA≒F)にする必要がある。

　　　ベクタ長が短いとs/cの精度が悪いのでaの精度も悪くなる。

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

#include <ap_int.h>

#include <cmath>  // atan()

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

 * atan2(s,c)

*  param range: FA<=30. F should be comparable to FA.

 *  input range: -2**(I-2)<=c,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;

  double e_min_c_val,e_min_s_val;

  double 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;

/* cos,sinをPI/M刻みで振ってaの最大誤差を表示する */

  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;

}

---

[←]　[↑]　[→]