Main Content

固定小数点の逆正接計算

スクリプトを開く

この例では、CORDIC アルゴリズム、多項式近似およびルックアップ テーブルというアプローチをそれぞれ使用して、固定小数点の 4 象限逆正接を計算する方法を説明します。これらの実装は、MATLAB® の組み込み関数 atan2 を近似しています。角度を推定する効率的な固定小数点の逆正接アルゴリズムは、ロボット工学の制御や無線通信における周波数追跡などの多数のアプリケーションで重要です。

CORDIC アルゴリズムを使用した atan2(y,x) の計算

はじめに

関数 cordicatan2 は CORDIC ベースのアルゴリズムを使用して MATLAB® 関数 atan2 を近似します。CORDIC は、COordinate Rotation DIgital Computer の略語です。ギブンス回転に基づく CORDIC アルゴリズム ([1] および [2] を参照) は、反復的なシフト加算演算だけが必要であるため、ハードウェア効率が最も高いアルゴリズムのうちの 1 つです。CORDIC アルゴリズムは、明示的な乗数を必要とせず、正弦関数、余弦関数、逆正弦関数、逆余弦関数、逆正接関数、ベクトル振幅関数、除算関数、平方根関数、双曲線関数、対数関数など、さまざまな関数の計算に適しています。

CORDIC のベクトル計算モード

CORDIC のベクトル モードの方程式は atan(y/x) の計算に広く使用されています。ベクトル モードでは、CORDIC 回転子が入力ベクトルを正の X 軸に対して回転させて、残差ベクトルの $$ y $$ 成分を最小化します。残差ベクトルの $$ y $$ 座標が正である場合、CORDIC 回転子は反復ごとに時計回りに (負の角度を使って) 回転し、それ以外の場合は反時計回りに (正の角度を使って) 回転します。反復の終了時に、角度のアキュムレータが 0 に初期化された場合、累計された回転角は元の入力ベクトルの角度です。

ベクトル モードでは、CORDIC 方程式は次のようになります。

$$ x_{i+1} = x_{i} - y_{i}*d_{i}*2^{-i} $$

$$ y_{i+1} = y_{i} + x_{i}*d_{i}*2^{-i} $$

$$ z_{i+1} = z_{i} + d_{i}*atan(2^{-i}) $$ は角度のアキュムレータです。

ここで、$$ y_{i} < 0 $$ の場合は $$ d_{i} = +1 $$、それ以外の場合は $$ -1 $$ です。

$$ i = 0, 1, ..., N-1 $$ および $$ N $$ は、反復の総数です。

$$ N $$ が $$ +\infty $$ に近づくにつれて、次のようになります。

$$ x_{N} = A_{N}\sqrt{x_{0}^2+y_{0}^2} $$

$$ y_{N} = 0 $$

$$ z_{N} = z_{0} + atan(y_{0}/x_{0}) $$

$$ A_{N} =&#xA;1/(cos(atan(2^{0}))*cos(atan(2^{-1}))*...*cos(atan(2^{-(N-1)})))&#xA; = \prod_{i=0}^{N-1}{\sqrt{1+2^{-2i}}}&#xA; $$

上記で説明しているように、逆正接は、角度のアキュムレータが 0 に初期化されているベクトル モードの CORDIC 回転子、つまり $$ z_{0}=0, $$ と $$ z_{N} \approx atan(y_{0}/x_{0}) $$ を使用して直接計算できます。

CORDICATAN2 のコードについて

はじめに

関数 cordicatan2 は x および y の要素の 4 象限逆正接を計算します ($$ -\pi \leq ATAN2(y,x) \leq +\pi $$)。cordicatan2 はベクトル モードの CORDIC アルゴリズムを使用し、上記の CORDIC 方程式に従って逆正接を計算します。

初期化

関数 cordicatan2 は以下の初期化手順を実行します。

  • $$ x_{0} $$ を X の初期入力値に設定します。

  • $$ y_{0} $$ を Y の初期入力値に設定します。

  • $$ z_{0} $$ をゼロに設定します。

$$ N $$ 回の反復の後、これらの初期値によって $$ z_{N} \approx atan(y_{0}/x_{0}) $$ が導かれます。

固定小数点および浮動小数点の共有 CORDIC カーネル コード

CORDIC アルゴリズム (ベクトル モード) のカーネル部分の MATLAB コードは以下のとおりです (スカラー x、y および z の場合)。固定小数点と浮動小数点の両方の演算には同じコードが使用されます。

function [x, y, z] = cordic_vectoring_kernel(x, y, z, inpLUT, n)
% Perform CORDIC vectoring kernel algorithm for N kernel iterations.
xtmp = x;
ytmp = y;
for idx = 1:n
    if y < 0
        x(:) = x - ytmp;
        y(:) = y + xtmp;
        z(:) = z - inpLUT(idx);
    else
        x(:) = x + ytmp;
        y(:) = y - xtmp;
        z(:) = z + inpLUT(idx);
    end
    xtmp = bitsra(x, idx); % bit-shift-right for multiply by 2^(-idx)
    ytmp = bitsra(y, idx); % bit-shift-right for multiply by 2^(-idx)
end

ベクトル モードの CORDIC の反復の可視化

CORDIC アルゴリズムは、通常、指定された回数 (定数) の反復を実行します。これは、CORDIC 反復を早期に終了するとパイプライン化されたコードが分断されて $$ n $$ が変化するため、CORDIC ゲイン $$ A_{n} $$ が一定でなくなるためです。

$$ n $$ が非常に大きな値の場合、CORDIC アルゴリズムが収束することは保証されますが、常に単調であるとは限りません。次の例からわかるように、中間の反復がベクトルを次の反復よりも正の X 軸に近づけるように回転する場合があります。通常、総反復回数を増やすと精度を上げることができます。

例

次の例では、反復 5 は反復 6 より優れた角度の推定を提供し、CORDIC アルゴリズムは後半の反復で収束します。

角度が $$ \theta = 43 $$ 度で、振幅が 1 である入力ベクトルを初期化します。

origFormat = get(0, 'format'); % store original format setting;
                               % restore this setting at the end.
format short
%
theta = 43*pi/180;  % input angle in radians
Niter = 10;         % number of iterations
inX   = cos(theta); % x coordinate of the input vector
inY   = sin(theta); % y coordinate of the input vector
%
% pre-allocate memories
zf = zeros(1, Niter);
xf = [inX, zeros(1, Niter)];
yf = [inY, zeros(1, Niter)];
angleLUT = atan(2.^-(0:Niter-1)); % pre-calculate the angle lookup table
%
% Call CORDIC vectoring kernel algorithm
for k = 1:Niter
   [xf(k+1), yf(k+1), zf(k)] = fixed.internal.cordic_vectoring_kernel_private(inX, inY, 0, angleLUT, k);
end

次の出力は、反復を 10 回行った CORDIC 角度の累積 (単位は度) を示します。5 回目の反復で生じた誤差は 6 回目の反復より少なく、計算された角度がその後の実際の入力角度に高速で収束したことに注意してください。

angleAccumulator = zf*180/pi; angleError = angleAccumulator - theta*180/pi;
fprintf('Iteration: %2d, Calculated angle: %7.3f, Error in degrees: %10g, Error in bits: %g\n',...
        [(1:Niter); angleAccumulator(:)'; angleError(:)';log2(abs(zf(:)'-theta))]);

Iteration:  1, Calculated angle:  45.000, Error in degrees:          2, Error in bits: -4.84036
Iteration:  2, Calculated angle:  18.435, Error in degrees:   -24.5651, Error in bits: -1.22182
Iteration:  3, Calculated angle:  32.471, Error in degrees:   -10.5288, Error in bits: -2.44409
Iteration:  4, Calculated angle:  39.596, Error in degrees:   -3.40379, Error in bits: -4.07321
Iteration:  5, Calculated angle:  43.173, Error in degrees:   0.172543, Error in bits: -8.37533
Iteration:  6, Calculated angle:  41.383, Error in degrees:   -1.61737, Error in bits: -5.14671
Iteration:  7, Calculated angle:  42.278, Error in degrees:  -0.722194, Error in bits: -6.3099
Iteration:  8, Calculated angle:  42.725, Error in degrees:   -0.27458, Error in bits: -7.70506
Iteration:  9, Calculated angle:  42.949, Error in degrees: -0.0507692, Error in bits: -10.1403
Iteration: 10, Calculated angle:  43.061, Error in degrees:  0.0611365, Error in bits: -9.87218

N が $$ +\infty $$ に近づくにつれて、CORDIC 回転子のゲイン $$ A_{N} $$ は 1.64676 に近づきます。この例では、入力 $$ (x_{0},y_{0}) $$ は単位円にあったため、回転子の初期振幅は 1 です。次の出力は、10 回の反復における回転子振幅を示します。

rotatorMagnitude = sqrt(xf.^2+yf.^2); % CORDIC rotator gain through iterations
fprintf('Iteration: %2d, Rotator magnitude: %g\n',...
    [(0:Niter); rotatorMagnitude(:)']);

Iteration:  0, Rotator magnitude: 1
Iteration:  1, Rotator magnitude: 1.41421
Iteration:  2, Rotator magnitude: 1.58114
Iteration:  3, Rotator magnitude: 1.6298
Iteration:  4, Rotator magnitude: 1.64248
Iteration:  5, Rotator magnitude: 1.64569
Iteration:  6, Rotator magnitude: 1.64649
Iteration:  7, Rotator magnitude: 1.64669
Iteration:  8, Rotator magnitude: 1.64674
Iteration:  9, Rotator magnitude: 1.64676
Iteration: 10, Rotator magnitude: 1.64676

$$ \sqrt{x_{0}^{2} + y_{0}^{2}} = 1 $$ であるため、$y_{n}$ が 0 に近づくと $x_{n}$ は $$ A_{n} \sqrt{x_{0}^{2} + y_{0}^{2}} = A_{n}, $$ に近づきます。

y_n = yf(end)

y_n =

   -0.0018


x_n = xf(end)

x_n =

    1.6468


figno = 1;
fidemo.fixpt_atan2_demo_plot(figno, xf, yf) %Vectoring Mode CORDIC Iterations


figno = figno + 1; %Cumulative Angle and Rotator Magnitude Through Iterations
fidemo.fixpt_atan2_demo_plot(figno,Niter, theta, angleAccumulator, rotatorMagnitude)

CORDIC アルゴリズムの全体的な誤差解析の実行

全体的な誤差は、次の 2 つの部分から構成されます。

  1. アルゴリズム誤差。基本的な角度の有限数で表されている CORDIC 回転角度から生じます。

  2. 量子化誤差または丸め誤差。角度ルックアップ テーブルの有限の精度表現や固定小数点演算で使用される有限の精度演算から生じます。

CORDIC アルゴリズム誤差の計算

theta  = (-178:2:180)*pi/180; % angle in radians
inXflt = cos(theta); % generates input vector
inYflt = sin(theta);
Niter  = 12; % total number of iterations
zflt   = cordicatan2(inYflt, inXflt, Niter); % floating-point results

CORDIC 計算を組み込み関数 atan2 と比較することによって、CORDIC アルゴリズム誤差の最大値を計算します。

format long
cordic_algErr_real_world_value = max(abs((atan2(inYflt, inXflt) - zflt)))

cordic_algErr_real_world_value =

     4.753112306290497e-04

対数の底が 2 の誤差は、反復回数に関連しています。この例では、反復を 12 回行った (つまり、11 個の 2 進数の精度になった) ため、誤差の大きさは $$ 2^{-11} $$ よりも小さくなります。

cordic_algErr_bits = log2(cordic_algErr_real_world_value)

cordic_algErr_bits =

 -11.038839889583048

"反復回数と精度の関係"

量子化誤差が全体的な誤差の大半を占める場合、つまり、量子化誤差がアルゴリズム誤差よりも多いときは、反復の合計回数を増加しても固定小数点の CORDIC アルゴリズムの全体的な誤差はあまり減りません。量子化誤差がアルゴリズム誤差よりも少なくなるように、小数部の長さと反復の合計回数を選択しなければなりません。CORDIC アルゴリズムでは、精度は反復ごとに 1 ビットずつ上がります。したがって、入力データの精度を超える多数の反復を選択する理由はありません。

反復回数と精度の関係は、アルゴリズムの右シフト ステップでもわかります。以下に、上記の反時計回りの回転の例を示します。

x(:) = x0 - bitsra(y,i);
y(:) = y + bitsra(x0,i);

i が y および x0 の語長と等しい場合、bitsra(y,i) と bitsra(x0,i) は 0 までシフトし、次のステップに何ら関与しません。

固定小数点アルゴリズムからの誤差を測定し、入力値の差分を測定しないように、浮動小数点参照を固定小数点 CORDIC アルゴリズムと同じ入力で計算します。

inXfix = sfi(inXflt, 16, 14);
inYfix = sfi(inYflt, 16, 14);
zref   = atan2(double(inYfix), double(inXfix));
zfix8  = cordicatan2(inYfix, inXfix, 8);
zfix10 = cordicatan2(inYfix, inXfix, 10);
zfix12 = cordicatan2(inYfix, inXfix, 12);
zfix14 = cordicatan2(inYfix, inXfix, 14);
zfix15 = cordicatan2(inYfix, inXfix, 15);
cordic_err = bsxfun(@minus,zref,double([zfix8;zfix10;zfix12;zfix14;zfix15]));

誤差は、反復回数と入力データの精度によって異なります。上記の例では、入力データは [-1, +1] の範囲にあり、小数部の長さは 14 です。各反復における最大誤差を示している次の表と、CORDIC アルゴリズムの全体的な誤差を示している図を見ると、データの精度に達するまで、誤差が反復ごとに約 1 ビットずつ減少しているのがわかります。

iterations = [8, 10, 12, 14, 15];
max_cordicErr_real_world_value = max(abs(cordic_err'));
fprintf('Iterations: %2d, Max error in real-world-value: %g\n',...
    [iterations; max_cordicErr_real_world_value]);

Iterations:  8, Max error in real-world-value: 0.00773633
Iterations: 10, Max error in real-world-value: 0.00187695
Iterations: 12, Max error in real-world-value: 0.000501175
Iterations: 14, Max error in real-world-value: 0.000244621
Iterations: 15, Max error in real-world-value: 0.000244621

max_cordicErr_bits = log2(max_cordicErr_real_world_value);
fprintf('Iterations: %2d, Max error in bits: %g\n',[iterations; max_cordicErr_bits]);

Iterations:  8, Max error in bits: -7.01414
Iterations: 10, Max error in bits: -9.05739
Iterations: 12, Max error in bits: -10.9624
Iterations: 14, Max error in bits: -11.9972
Iterations: 15, Max error in bits: -11.9972

figno = figno + 1;
fidemo.fixpt_atan2_demo_plot(figno, theta, cordic_err)

FIACCEL を使用した固定小数点アルゴリズム CORDICATAN2 の高速化

MATLAB® コマンド fiaccel を使用して、MATLAB コードから MEX 関数を生成できます。通常は、生成した MEX 関数を実行するとシミュレーションの速度が向上しますが、実際にどの程度速度が向上するかは、使用されるシミュレーション プラットフォームによります。以下の例は、fiaccel を使用して固定小数点アルゴリズム cordicatan2 を高速化する方法を示します。

関数 fiaccel は、MATLAB コードを MEX 関数にコンパイルします。このステップには、一時ディレクトリの作成とそのディレクトリにおける書き込み権限が必要です。

tempdirObj = fidemo.fiTempdir('fixpt_atan2_demo');

coder.newtype('constant',12) を使用して反復回数が定数 (たとえば、12) になるように宣言すると、コンパイルされた角度ルックアップ テーブルも定数になるため、各反復で計算されません。また、コンパイルされた MEX ファイル cordicatan2_mex を呼び出すときには、反復回数に対する入力引数を指定する必要がありません。反復回数を渡すと、MEX 関数はエラーになります。

入力パラメーターのデータ型は、関数 cordicatan2 で固定小数点または浮動小数点の計算を行うかどうかを決定します。MATLAB がこのファイルに対してコードを生成すると、コードは特定のデータ型に対してのみ生成されます。たとえば、入力が固定小数点である場合は、固定小数点コードのみが生成されます。

inp = {inYfix, inXfix, coder.newtype('constant',12)}; % example inputs for the function
fiaccel('cordicatan2', '-o', 'cordicatan2_mex',  '-args', inp)

最初に、cordicatan2 を呼び出して 4 象限 atan2 のベクトルを計算します。

tstart = tic;
cordicatan2(inYfix,inXfix,Niter);
telapsed_Mcordicatan2 = toc(tstart);

次に、MEX 関数 cordicatan2_mex を呼び出して、4 象限 atan2 のベクトルを計算します。

cordicatan2_mex(inYfix,inXfix); % load the MEX file
tstart = tic;
cordicatan2_mex(inYfix,inXfix);
telapsed_MEXcordicatan2 = toc(tstart);

ここで、速度を比較します。MATLAB コマンド ウィンドウに次のコマンドを入力して、特定のプラットフォームにおける速度の改善を確認します。

fiaccel_speedup = telapsed_Mcordicatan2/telapsed_MEXcordicatan2;

一時ディレクトリをクリーンアップするには、次のコマンドを実行します。

clear cordicatan2_mex;
status = tempdirObj.cleanUp;

チェビシェフ多項式近似を使用した atan2(y,x) の計算

多項式近似は、乗算累積 (MAC) 中心のアルゴリズムです。atan(x) のような非線形関数の DSP 実装に適しています。

特定の次元の多項式と [-1, +1] の区間で評価される特定の関数 f(x) = atan(x) に対して、多項式近似理論は $$ |P(x)-f(x)| $$ の最大値を最小化する多項式を求めようとします。ここで、P(x) は近似多項式です。一般的に、チェビシェフ多項式で与えられた関数を近似し、多項式を目的の次元でカット オフすることによって、最適な多項式に非常に近い多項式を得ることができます。

第 1 種チェビシェフ多項式を使用した区間 [-1, +1] における逆正接の近似は、次の式にまとめられます。

$$ atan(x) = 2\sum_{n=0}^{\infty} {(-1)^{n}q^{2n+1} \over (2n+1)}&#xA;T_{2n+1}(x) $$

ここで、

$$ q = 1/(1+\sqrt{2}) $$

$$ x \in [-1, +1] $$

$$ T_{0}(x) = 1 $$

$$ T_{1}(x) = x $$

$$ T_{n+1}(x) = 2xT_{n}(x) - T_{n-1}(x). $$

したがって、3 次のチェビシェフ多項式近似は次のとおりです。

$$ atan(x) = 0.970562748477141*x - 0.189514164974601*x^{3}. $$

5 次のチェビシェフ多項式近似は次のとおりです。

$$ atan(x) = 0.994949366116654*x - 0.287060635532652*x^{3}&#xA; + 0.078037176446441*x^{5}. $$

7 次のチェビシェフ多項式近似は次のとおりです。

$$ \begin{array}{lllll}&#xA; atan(x) &#38; = &#38; 0.999133448222780*x &#38; - &#38; 0.320533292381664*x^{3} \\&#xA; &#38; + &#38; 0.144982490144465*x^{5} &#38; - &#38; 0.038254464970299*x^{7}.&#xA;\end{array} $$

逆正接関数のプロパティに基づいた角度補正を使用して 4 象限出力を得ることができます。

CORDIC アルゴリズムと多項式近似アルゴリズムのアルゴリズム誤差の比較

一般的に、多項式近似の次数が高いほどより正確な最終結果が得られます。ただし、多項式近似の次数が高いとアルゴリズムの複雑さも増すため、MAC 演算数が増えて、より多くのメモリが必要になります。CORDIC アルゴリズムおよび MATLAB の関数 atan2 と矛盾しないようにするために、入力引数は y/x の比ではなく x 座標と y 座標の両方で構成されます。

量子化誤差をなくすには、CORDIC およびチェビシェフ多項式近似アルゴリズムの浮動小数点の実装を以下の比較で使用します。アルゴリズム誤差の比較から、CORDIC の反復回数を増やすと誤差が減ることがわかります。また、反復を 12 回行った CORDIC アルゴリズムでは、5 次のチェビシェフ多項式近似より角度推定値がわずかに優れていることもわかります。3 次のチェビシェフ多項式の近似誤差は、5 次のチェビシェフ多項式の 8 倍の大きさです。多項式の順序や次数は、角度推定値の必要な精度とハードウェア制約に基づいて選択しなければなりません。

atan(x) に対するチェビシェフ多項式近似の係数は、昇順の x で示されます。

constA3 = [0.970562748477141, -0.189514164974601]; % 3rd order
constA5 = [0.994949366116654,-0.287060635532652,0.078037176446441]; % 5th order
constA7 = [0.999133448222780 -0.320533292381664 0.144982490144465...
          -0.038254464970299]; % 7th order

theta   = (-90:1:90)*pi/180; % angle in radians
inXflt  = cos(theta);
inYflt  = sin(theta);
zfltRef = atan2(inYflt, inXflt); % Ideal output from ATAN2 function
zfltp3  = fidemo.poly_atan2(inYflt,inXflt,3,constA3); % 3rd order polynomial
zfltp5  = fidemo.poly_atan2(inYflt,inXflt,5,constA5); % 5th order polynomial
zfltp7  = fidemo.poly_atan2(inYflt,inXflt,7,constA7); % 7th order polynomial
zflt8   = cordicatan2(inYflt, inXflt,  8); % CORDIC alg with 8 iterations
zflt12  = cordicatan2(inYflt, inXflt, 12); % CORDIC alg with 12 iterations

以下に、反復を 8 回と 12 回行った CORDIC アルゴリズムのアルゴリズム誤差の最大値 (またはアルゴリズム誤差の無限大ノルム) を示します。

cordic_algErr    = [zfltRef;zfltRef] - [zflt8;zflt12];
max_cordicAlgErr = max(abs(cordic_algErr'));
fprintf('Iterations: %2d, CORDIC algorithmic error in real-world-value: %g\n',...
    [[8,12]; max_cordicAlgErr(:)']);

Iterations:  8, CORDIC algorithmic error in real-world-value: 0.00772146
Iterations: 12, CORDIC algorithmic error in real-world-value: 0.000483258

対数の底が 2 の誤差は、精度の 2 進数の数を示します。CORDIC アルゴリズムの 12 回目の反復で推定角度精度は $$ 2^{-11} $$ です。

max_cordicAlgErr_bits = log2(max_cordicAlgErr);
fprintf('Iterations: %2d, CORDIC algorithmic error in bits: %g\n',...
    [[8,12]; max_cordicAlgErr_bits(:)']);

Iterations:  8, CORDIC algorithmic error in bits: -7.01691
Iterations: 12, CORDIC algorithmic error in bits: -11.0149

次のコードは、3 次、5 次、および 7 次の多項式近似のアルゴリズム誤差の最大値を示します。

poly_algErr    = [zfltRef;zfltRef;zfltRef] - [zfltp3;zfltp5;zfltp7];
max_polyAlgErr = max(abs(poly_algErr'));
fprintf('Order: %d, Polynomial approximation algorithmic error in real-world-value: %g\n',...
    [3:2:7; max_polyAlgErr(:)']);

Order: 3, Polynomial approximation algorithmic error in real-world-value: 0.00541647
Order: 5, Polynomial approximation algorithmic error in real-world-value: 0.000679384
Order: 7, Polynomial approximation algorithmic error in real-world-value: 9.16204e-05

対数の底が 2 の誤差は、精度の 2 進数の数を示します。

max_polyAlgErr_bits = log2(max_polyAlgErr);
fprintf('Order: %d, Polynomial approximation algorithmic error in bits: %g\n',...
    [3:2:7; max_polyAlgErr_bits(:)']);

Order: 3, Polynomial approximation algorithmic error in bits: -7.52843
Order: 5, Polynomial approximation algorithmic error in bits: -10.5235
Order: 7, Polynomial approximation algorithmic error in bits: -13.414

figno = figno + 1;
fidemo.fixpt_atan2_demo_plot(figno, theta, cordic_algErr, poly_algErr)

浮動小数点チェビシェフ多項式近似アルゴリズムの固定小数点への変換

入力および出力の語長がハードウェアによって 16 ビットに制限されており、5 次のチェビシェフ多項式を近似で使用すると仮定します。入力 x、y、および y/x のダイナミック レンジはすべて [-1, +1] 内であるため、語長が 16 ビットで、小数部の長さが 14 ビットである符号付き固定小数点入力データ型を選択するとオーバーフローを回避できます。多項式の係数は純粋に分数で (-1, +1) 内であるため、語長が 16 ビットで、小数部の長さが 15 ビット (最高の精度) である符号付き固定小数点として、それらのデータ型を選択できます。$$ (y/x)^{n} $$ が [-1, +1] 内であり、係数の乗算および $$ (y/x)^{n} $$ が (-1, +1) 内であるため、アルゴリズムはロバストです。したがって、ダイナミック レンジは拡大せず、事前定義した固定小数点データ型により、オーバーフローは発生しません。

CORDIC アルゴリズムと同様に、4 象限多項式近似ベースの atan2 アルゴリズムは $$ [-\pi, \pi] $$ 内の推定角度を出力します。したがって、13 ビットの出力の小数部の長さを選択するとオーバーフローを回避して [-4, +3.9998779296875] のダイナミック レンジを指定できます。

区間 [-1, +1] における逆正接の基本的な浮動小数点チェビシェフ多項式近似はファイル poly_atan2.m 内でローカル関数 chebyPoly_atan_fltpt として実装されています。

   function z = chebyPoly_atan_fltpt(y,x,N,constA,Tz,RoundingMethodStr)
   tmp = y/x;
   switch N
       case 3
           z = constA(1)*tmp + constA(2)*tmp^3;
       case 5
           z = constA(1)*tmp + constA(2)*tmp^3 + constA(3)*tmp^5;
       case 7
           z = constA(1)*tmp + constA(2)*tmp^3 + constA(3)*tmp^5 + constA(4)*tmp^7;
       otherwise
           disp('Supported order of Chebyshev polynomials are 3, 5 and 7');
   end

区間 [-1, +1] における逆正接の基本的な固定小数点チェビシェフ多項式近似はファイル poly_atan2.m 内でローカル関数 chebyPoly_atan_fixpt として実装されています。

   function z = chebyPoly_atan_fixpt(y,x,N,constA,Tz,RoundingMethodStr)
   z = fi(0,'numerictype', Tz, 'RoundingMethod', RoundingMethodStr);
   Tx = numerictype(x);
   tmp = fi(0, 'numerictype',Tx, 'RoundingMethod', RoundingMethodStr);
   tmp(:) = Tx.divide(y, x); % y/x;
   tmp2 = fi(0, 'numerictype',Tx, 'RoundingMethod', RoundingMethodStr);
   tmp3 = fi(0, 'numerictype',Tx, 'RoundingMethod', RoundingMethodStr);
   tmp2(:) = tmp*tmp;  % (y/x)^2
   tmp3(:) = tmp2*tmp; % (y/x)^3
   z(:) = constA(1)*tmp + constA(2)*tmp3; % for order N = 3
   if (N == 5) || (N == 7)
       tmp5 = fi(0, 'numerictype',Tx, 'RoundingMethod', RoundingMethodStr);
       tmp5(:) = tmp3 * tmp2; % (y/x)^5
       z(:) = z + constA(3)*tmp5; % for order N = 5
       if N == 7
           tmp7 = fi(0, 'numerictype',Tx, 'RoundingMethod', RoundingMethodStr);
           tmp7(:) = tmp5 * tmp2; % (y/x)^7
           z(:) = z + constA(4)*tmp7; %for order N = 7
       end
   end

チェビシェフ多項式近似を使用した普遍的な 4 象限 atan2 計算は、poly_atan2.m ファイルに実装されます。

   function z = poly_atan2(y,x,N,constA,Tz,RoundingMethodStr)
   if nargin < 5
       % floating-point algorithm
       fhandle = @chebyPoly_atan_fltpt;
       Tz = [];
       RoundingMethodStr = [];
       z = zeros(size(y));
   else
       % fixed-point algorithm
       fhandle = @chebyPoly_atan_fixpt;
       %pre-allocate output
       z = fi(zeros(size(y)), 'numerictype', Tz, 'RoundingMethod', RoundingMethodStr);
   end
   % Apply angle correction to obtain four quadrant output
   for idx = 1:length(y)
      % first quadrant
      if abs(x(idx)) >= abs(y(idx))
          % (0, pi/4]
          z(idx) = feval(fhandle, abs(y(idx)), abs(x(idx)), N, constA, Tz, RoundingMethodStr);
      else
          % (pi/4, pi/2)
          z(idx) = pi/2 - feval(fhandle, abs(x(idx)), abs(y(idx)), N, constA, Tz, RoundingMethodStr);
      end
      if x(idx) < 0
          % second and third quadrant
          if y(idx) < 0
              z(idx) = -pi + z(idx);
          else
             z(idx) = pi - z(idx);
          end
      else % fourth quadrant
          if y(idx) < 0
              z(idx) = -z(idx);
          end
      end
   end

多項式近似アルゴリズムの全体的な誤差解析の実行

CORDIC アルゴリズムと同様に、多項式近似アルゴリズムの全体的な誤差は 2 つの部分で構成されます。すなわち、アルゴリズム誤差と量子化誤差です。多項式近似のアルゴリズム誤差は、前の節で解析し、CORDIC アルゴリズムのアルゴリズム誤差と比較しました。

量子化誤差の計算

量子化誤差を、固定小数点多項式近似を浮動小数点多項式近似と比較することで計算します。

最も近い偶数方向への丸めを使用して入力と係数を量子化します。

inXfix = fi(fi(inXflt,  1, 16, 14,'RoundingMethod','Convergent'),'fimath',[]);
inYfix = fi(fi(inYflt,  1, 16, 14,'RoundingMethod','Convergent'),'fimath',[]);
constAfix3 = fi(fi(constA3, 1, 16,'RoundingMethod','Convergent'),'fimath',[]);
constAfix5 = fi(fi(constA5, 1, 16,'RoundingMethod','Convergent'),'fimath',[]);
constAfix7 = fi(fi(constA7, 1, 16,'RoundingMethod','Convergent'),'fimath',[]);

Floor 丸めを使用して量子化誤差の最大値を計算します。

ord    = 3:2:7; % using 3rd, 5th, 7th order polynomials
Tz     = numerictype(1, 16, 13); % output data type
zfix3p = fidemo.poly_atan2(inYfix,inXfix,ord(1),constAfix3,Tz,'Floor'); % 3rd order
zfix5p = fidemo.poly_atan2(inYfix,inXfix,ord(2),constAfix5,Tz,'Floor'); % 5th order
zfix7p = fidemo.poly_atan2(inYfix,inXfix,ord(3),constAfix7,Tz,'Floor'); % 7th order
poly_quantErr = bsxfun(@minus, [zfltp3;zfltp5;zfltp7], double([zfix3p;zfix5p;zfix7p]));
max_polyQuantErr_real_world_value = max(abs(poly_quantErr'));
max_polyQuantErr_bits = log2(max_polyQuantErr_real_world_value);
fprintf('PolyOrder: %2d, Quant error in bits: %g\n',...
    [ord; max_polyQuantErr_bits]);

PolyOrder:  3, Quant error in bits: -12.7101
PolyOrder:  5, Quant error in bits: -12.325
PolyOrder:  7, Quant error in bits: -11.8416

全体的な誤差の計算

固定小数点多項式近似を組み込み関数 atan2 と比較することで、全体的な誤差を計算します。理想の参照出力は zfltRef です。7 次の多項式近似の全体的な誤差は、ほとんどが量子化誤差であり、固定小数点の算術演算からの入力データ、係数および丸め効果の有限精度が原因となります。

poly_err = bsxfun(@minus, zfltRef, double([zfix3p;zfix5p;zfix7p]));
max_polyErr_real_world_value = max(abs(poly_err'));
max_polyErr_bits = log2(max_polyErr_real_world_value);
fprintf('PolyOrder: %2d, Overall error in bits: %g\n',...
    [ord; max_polyErr_bits]);

PolyOrder:  3, Overall error in bits: -7.51907
PolyOrder:  5, Overall error in bits: -10.2497
PolyOrder:  7, Overall error in bits: -11.5883

figno = figno + 1;
fidemo.fixpt_atan2_demo_plot(figno, theta, poly_err)

"多項式近似における丸めモードの影響"

角度のアキュムレータで、反復 12 回、小数部の長さ 13 ビットの CORDIC アルゴリズムと比較すると、5 次のチェビシェフ多項式近似からは同様の順序の量子化誤差が得られます。次の例では、Nearest、Round、および Convergent 丸めモードの量子化誤差は Floor 丸めモードよりも小さくなります。

Floor 丸めを使用した量子化誤差の最大値は次のようになります。

poly5_quantErrFloor = max(abs(poly_quantErr(2,:)));
poly5_quantErrFloor_bits = log2(poly5_quantErrFloor)

poly5_quantErrFloor_bits =

 -12.324996933210334

比較する場合は、次のように Nearest 丸めを使用して量子化誤差の最大値を計算します。

zfixp5n = fidemo.poly_atan2(inYfix,inXfix,5,constAfix5,Tz,'Nearest');
poly5_quantErrNearest = max(abs(zfltp5 - double(zfixp5n)));
poly5_quantErrNearest_bits = log2(poly5_quantErrNearest)
set(0, 'format', origFormat); % reset MATLAB output format

poly5_quantErrNearest_bits =

 -13.175966487895451

ルックアップ テーブルを使用した atan2(y,x) の計算

固定小数点の逆正接近似を実装するために使用できるルックアップ テーブルに基づくアプローチは数多く存在します。以下に示すのは、単一の実数値ルックアップ テーブルと単純な最近傍線形内挿に基づく低コストのアプローチです。

単一のルックアップ テーブルに基づくアプローチ

Fixed-Point Designer™ に含まれる fi オブジェクトの atan2 メソッドは、単一のルックアップ テーブルに基づくアプローチを値間の単純な最近傍線形内挿と共に使用して、MATLAB® の組み込み浮動小数点関数 atan2 を近似します。このアプローチでは、実数値ルックアップ テーブルを小規模にすることができ、使用する算術も簡単です。

単一の実数値ルックアップ テーブルを使用することで、インデックスの計算と、非常に高い精度の結果を得るために必要な計算全体が単純化されます。これらの単純化により、比較的高速のパフォーマンスが実現するとともに、必要なメモリも比較的少なくなります。

ルックアップ テーブルに基づく ATAN2 の実装について

ルックアップ テーブルのサイズと精度

ルックアップ テーブルの設計で考慮するべき 2 つの重要な点は、サイズと精度です。考え得るすべての入力値 $$ y/x $$ に対応できる table を作成するのは不可能です。また、ルックアップ テーブル値の量子化のため、完璧に正確であることも不可能です。

妥協案として、Fixed-Point Designer の fi オブジェクトの atan2 メソッドは、実装の一部として 8 ビットのルックアップ テーブルを使用します。8 ビット テーブルの長さはわずか 256 要素なので、小さくて効率的です。さらに、多くのプラットフォームにおける 8 ビットは 1 バイトまたは 1 ワードのサイズと同じです。線形内挿と組み合わせて使用し、出力 (ルックアップ テーブル値) の精度を 16 ビットにすることで、8 ビットでアドレス指定可能なルックアップ テーブルは非常に高い精度とパフォーマンスを提供します。

アルゴリズム実装の概要

Fixed-Point Designer の実装をより深く理解するために、最初に 4 象限関数 atan2(y,x) の対称性を考慮します。常に x-y 空間の最初の八分円で (つまり、0 ～ pi/4 ラジアンの角度で) 逆正接を計算する場合は、任意の y 値 および x 値の結果の角度に対して八分円補正を実行できます。

前処理の一環として、y と x の符号と相対的な大きさが考慮され、除算が実行されます。y と x の符号と大きさに基づいて、y/x、x/y、-y/x、-x/y、-y/-x、-x/-y のうち 1 つだけが計算されます。y と x の符号と大きさに関する事前の知識に基づいて、非負の純粋な非整数であることが保証される符号なしの結果が計算されます。この値には、符号なしの 16 ビット小数固定小数点型が使用されます。

その後、純粋な非整数符号なし固定小数点結果の符号なし格納整数表現の最上位ビット (MSB) 8 個が、範囲 0 ～ pi/4 ラジアンの角度の値が収められた 8 ビット (長さ 256) ルックアップ テーブル値に対する直接のインデックスとして使用されます。テーブル ルックアップが 2 回実行されます。1 回は計算されたテーブル インデックス位置 lutValBelow で、もう 1 回は次のインデックス位置 lutValAbove で実行されます。

idxUint8MSBs = bitsliceget(idxUFIX16, 16, 9);
zeroBasedIdx = int16(idxUint8MSBs);
lutValBelow  = FI_ATAN_LUT(zeroBasedIdx + 1);
lutValAbove  = FI_ATAN_LUT(zeroBasedIdx + 2);

idxUFIX16 の残りの 8 個の最下位ビット (LSB) は、2 つのテーブル値の間の内挿に使用されます。LSB 値は、8 ビットの小数データ型 rFracNT を伴う正規化されたスケーリング ファクターとして扱われます。

rFracNT      = numerictype(0,8,8); % fractional remainder data type
idxFrac8LSBs = reinterpretcast(bitsliceget(idxUFIX16,8,1), rFracNT);
rFraction    = idxFrac8LSBs;

2 つのルックアップ テーブル値を剰余 (rFraction) 値とともに使用して、単純な最近傍線形内挿が実行されます。実数の乗算を使用して、2 点間の重み付き差分を決定します。その結果、内挿された固定小数点数の結果を得るための単純な計算 (1 つの積と 2 つの和) が導かれます。

temp = rFraction * (lutValAbove - lutValBelow);
rslt = lutValBelow + temp;

最後に、y と x の元の符号と相対的大きさに基づき、単純な八分円補正ロジックおよび算術を使用して、出力結果が形成されます。最初の八分円 [0, pi/4] の角度値の結果が定数を使って加算または減算されて、八分円補正された角度の出力が形成されます。

ATAN2 を使用した固定小数点の逆正接の計算

固定小数点または浮動小数点入力を使って、関数 atan2 を直接呼び出すことができます。固定小数点 atan2 の実装には、ルックアップ テーブルに基づくアルゴリズムが使用されます。

zFxpLUT = atan2(inYfix,inXfix);

全体的な誤差の計算

固定小数点ルックアップ テーブルに基づく近似を組み込み関数 atan2 と比較することで、全体的な誤差を計算できます。理想の参照出力は zfltRef です。

lut_err = bsxfun(@minus, zfltRef, double(zFxpLUT));
max_lutErr_real_world_value = max(abs(lut_err'));
max_lutErr_bits = log2(max_lutErr_real_world_value);
fprintf('Overall error in bits: %g\n', max_lutErr_bits);

Overall error in bits: -12.6743

figno = figno + 1;
fidemo.fixpt_atan2_demo_plot(figno, theta, lut_err)

固定小数点実装間での全体的な誤差の比較

上記で行ったように、固定小数点近似を組み込み関数 atan2 と比較することで、全体的な誤差を計算できます。理想の参照出力は zfltRef です。

zfixCDC15      = cordicatan2(inYfix, inXfix, 15);
cordic_15I_err = bsxfun(@minus, zfltRef, double(zfixCDC15));
poly_7p_err    = bsxfun(@minus, zfltRef, double(zfix7p));
figno = figno + 1;
fidemo.fixpt_atan2_demo_plot(figno, theta, cordic_15I_err, poly_7p_err, lut_err)

固定小数点近似アルゴリズムのコストの比較

固定小数点 CORDIC アルゴリズムには次の演算が必要です。

  • 反復あたり 1 回のテーブル ルックアップ

  • 反復あたり 2 回のシフト

  • 反復あたり 3 回の加算

N 次の固定小数点チェビシェフ多項式近似アルゴリズムには、次の演算が必要です。

  • 1 回の除算

  • 乗算 (N+1) 回

  • 加算 (N-1)/2 回

最近傍線形内挿を伴う単純化された単一のルックアップ テーブル アルゴリズムには、次の演算が必要です。

  • 1 回の除算

  • 2 回のテーブル ルックアップ

  • 1 回の乗算

  • 2 回の加算

実際のアプリケーションでは、固定小数点の逆正接を計算するためのアルゴリズムの選択は、通常、必要な精度、コスト、およびハードウェア制約によって異なります。

close all; % close all figure windows

参考文献

  1. Jack E. Volder, The CORDIC Trigonometric Computing Technique, IRE Transactions on Electronic Computers, Volume EC-8, September 1959, pp. 330-334.

  2. Ray Andraka, A survey of CORDIC algorithm for FPGA based computers, Proceedings of the 1998 ACM/SIGDA sixth international symposium on Field programmable gate arrays, Feb. 22-24, 1998, pp. 191-200.

%#ok<*NOPTS>